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Abstract 


The  purpose  of  this  study  was  to  determine  the  effects 
of  convection  heat  transfer,  target  vapor  condensation,  and 
radiation  heat  transfer  on  the  observed  enhanced  thermal 
coupling  of  a  pulsed  laser  to  an  Aluminum  target.  The  need 
for  this  study  is  to  develop  a  base  from  which  the  laser 
parameters  needed  to  most  efficiently  couple  laser  energy 
to  a  metal  target  may  be  predicted. 

An  attempt  was  made  to  analytically  duplicate  a  set  of 
experimental  data  in  which  enhanced  thermal  coupling  was 
observed.  Analytical  relations  were  developed  to  compute 
the  heat  flux,  into  the  surface  of  an  Aluminum  species,  due 
to  convection  heat  transfer,  target  vapor  condensation,  and 
plasma  reradiation.  From  these  relations,  the  most 
significant  coupling  mechanism  was  identified  to  be  vapor 
condensation.  Convection^  heat  transfer  was  found  to  have  a 
less  significant  effect  on  enhanced  coupling,  and  plasma 
reradiation  was  found  to  be  negligible.  The  calculations 
were  performed  by  a  three-dimensional  Lagrangian 
hydrocode.  The  results  of  this  work  are  somewhat 
inconclusive  due  to  the  extensive  computer  resources  used 
by  the  hydrocode;  however,  the  relations  presented  here  do 
seem  to  give  promising  results.  The  work  should  be 
continued  to  investigate  the  effects  of  laser  and  target 
parameters  on  enhanced  thermal  coupling. 


ENHANCED  THERMAL  COUPLING  BY  A  REPETITIVELY 

PULSED  LASER 

I .  Introduction 

Much  work  has  been  done  in  past  years,  primarily  by 
researchers  at  the  United  States  Air  Force  Weapons 
Laboratory,  in  studying  the  effects  of  high-power  lasers  on 
metallic  surfaces.  In  particular,  much  research  has  been 
directed  towards  studying  the  thermal  coupling  of  lasers  to 
metals.  The  amount  of  thermal  coupling  is  a  direct 
indication  of  how  efficiently  laser  energy  is  being 
delivered  to  the  target.*  In  most  practical  applications  of 
laser  energy,  it  is  sometimes  desirable  to  deliver  the  most 
energy  to  the  target  for  the  least  investment,  in  terms  of 
total  laser  energy;  thus,  it  is  desirable  to  maximize  the 
amount  of  "thermal"  coupling.  Experimental  data,  though  the 
amount  is  relatively  limited,  is  available  on  the  subject. 
What  is  needed,  however,  is  a  theoretical  explanation  of  how 
the  amount  of  thermal  coupling  may  be  increased.  Since  such 
an  explanation  is  not  currently  available,  a  natural 
"starting  point"  for  its  development  is  an  attempt  to 
analytically  duplicate  a  set  of  experimental  data.  Such  a 
model,  due  to  the  fact  that  it  is  derived  from  a  given  set 
of  experimental  data,  will  not  be  expected  to  yield  accurate 

*  The  words  target,  metal,  and  surface  will  be  used 
synonymously  in  this  document. 
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results  for  all  experimental  conditions;  however,  in 
theory,  it  can  be  modified  to  include  the  effect  of 
experimental  parameters.  Once  the  model  is  fully 
developed,  it  will  then  be  possible  to  predict  the  laser 
parameters  needed  to  maximize  the  amount  of  thermal 
coupling  to  a  particular  type  of  target. 

It  is  known  {13  that,  under  certain  conditions,  the 
thermal  coupling  of  a  pulsed  laser  (with  average  intensity 
I)  to  a  target  is  an  order  of  magnitude  greater  than  that 
of  a  continuous  wave  laser  (with  intensity  I).  For 
example.  Hall  {2}**  reports  that  the  coupling  of  a  10.6 
micron  continuous  wave  (CW)  CO2  laser  to  cold  Aluminum  is 
about  3%;  whereas,  the  coupling  of  a  pulsed  CO2  laser  of 
the  same  wavelength  is  about  17.5%  for  a  9.3  J  pulse.  The 
as  yet  unanswered  question  is  then,  "why  does  the 
enhancement  occur?"  It  is  generally  accepted  that  the 
overall  reason  for  the  enhancement  is  that  much  higher 
intensities  are  attainable  when  a  laser  is  operated  in  a 
pulsed  mode,  as  opposed  to  a  continuous  wave  mode. 
Furthermore,  it  has  been  observed  [33  that  a  prerequisite 
for  enhancement  is  the  ignition  of  an  air  breakdown  plasma 
called  -  laser  supported  detonation  wave  (LSD) .  It 
appears,  therefore,  that  the  phenomenon  of  enhanced 
thermal  coupling  is  due  to  the  presence  of  an  LSD. 

**  Hall's  experimental  data  is  the  background  for  this  re¬ 
search,  and  is  discussed  in  detail  in  Section  II. 


The  LSD  is  a  highly  absorptive  plasma  front,  considered 
to  be  optically  thick.  For  this  reason,  it  has  been 
acknowledged  £4,  5,  6}  that  the  enhancement  is  due  to  the 
efficient  transfer  of  energy  from  the  plasma  to  the  target 
by  convection  and  radiation  heat  transfer.  Approximate 
heat  transfer  calculations,  however,  have  been  made  and 
they  fail  to  match  experimental  data.  For  example,  Stamm, 
Nielsen,  and  Jackson  £7)  report  that  radiation  is  not  the 
cause  of  the  enhancement.  Also,  Jumper  and  Jackson  [81 
report  that  convection  heat  transfer  is  probably  not 
responsible  for  the  enhancement.  There  is,  then,  a  large 
discrepancy  between  experimental  data  and  approximate  heat 
transfer  calculations. 

When  a  laser  pulse  of  sufficient  intensity  strikes  a 
target,  boiling  and  vaporization  of  the  target  may  occur. 

If  a  LSD  is  ignited,  it  propagates  up  the  laser  beam  and, 
as  it  passes  by  a  region  of  target  vapor,  the  high 
pressures  associated  with  the  LSD  force  the  target  vapor 
back  onto  the  surface  where  it  flows  out  radially  and 
condenses.  Jumper  £93  has  suggested  that  this  condensation 
may  be  responsible  for  the  enhancement  and  he  has  shown 
£10}  that  sufficient  energy  resides  in  the  condensing  vapor 
to  account  for  observed  values  of  enhanced  coupling.  This 
theory  is  worthy  of  further  investigation  since  the  vapor, 
when  it  condenses,  gives  up  not  only  the  latent  heat  of 
vaporization  required  to  vaporize  it,  but  also  the  energy 


it  gains  while  in  the  vapor  state.  Therefore,  it  appears 
that  if  a  LSD  forces  target  vapor  onto  the  surface  and  if 
some  or  all  the  vapor  condenses,  the  final  energy  residing 
in  the  target  will  be  greater  and  thus  the  amount  of 
thermal  coupling  will  be  enhanced. 

The  foremost  object  of  this  research  effort  was  to 
perform  detailed  calculations  of  the  convection  heat 
transfer  taking  place  at  the  surface  of  an  aluminum 
species,  under  the  presence  of  a  LSD.  The  results  of  the 
heat  transfer  calculations  were  compared  to  Hall's 
experimental  data  [2],  in  order  that  the  role  of  convection 
heat  transfer  on  enhanced  coupling  could  be  determined. 

The  second  object  of  this  research  was  to  examine  the  role 
of  vapor  condensation  on  enhanced  coupling.  The  vapor 
condensation  calculations  were  performed  in  more  detail 
than  those  made  by  Jumper  and  Jackson  [10).  The  third  and 
final  object  of  this  research  was  to  predict  the  role  of 
plasma  radiation  on  enhanced  coupling.  It  was  projected 
that  these  calculations  would  identify  the  most  important 
cause  of  enhanced  thermal  coupling  for  Hall's  9.3  J  pulse. 
The  effect  of  beam  parameters  on  coupling  was  not  studied 
here;  however,  it  is  hoped  that  the  models  developed  here 
might  eventually  be  extended  to  include  these  parameters. 
Future  study  in  this  area  is  needed. 

The  calculations  needed  for  this  study  were  performed 
numerically  by  a  three-dimensional  Lagrangian  hydrocode. 


given  by  Nielsen  [11).  The  hydrocode  is  explained  in 
detail  in  Section  III  and  it  suffices  to  say  here  that  the 
hydrocode  solves  the  inviscid  "control  volume"  equations 


for  the  conservation  of  mass,  momentum,  and  energy, 
augmented  with  the  perfect  gas  equations  of  state  for  air 
and  target  vapor.  The  hydrocode  gives  the  following 
time-dependent  flow  properties  for  air,  through  which  a 
time-varying,  spatial-varying  laser  beam  is  propagating: 
pressure,  total  density,  vapor  density,  energy  density, 
components  of  momentum  (velocity)  in  all  spatial  directions 
of  interest,  and  temperature.  The  hydrocode  may  be 
modified  to  include  the  presence  of  a  target  and,  of 
course,  since  this  problem  is  a  laser-target  interaction 
problem,  the  presence  of  a  target  was  assumed. 

The  reader  is  informed  that  the  results  of  this  study 
are  somewhat  inconclusive,  as  the  amount  of  numerical  data 
for  analysis  was  limited.  This  was  due  to  the  fact  that 
the  hydrocode,  because  of  its  size  and  multi¬ 
dimensionality,  requires  massive  amounts  of  computer  time 
and  computer  resources;  whereas  the  amount  of  computer  time 
allocated  to  this  study  was  limited.  The  large  computer 
time  requirement  arises  due  to  the  fact  that  this  study  was 
one  of  boundary  layer  effects  and,  in  order  to  accurately 
predict  the  properties  at  the  edge  of  a  boundary  layer,  of 
unknown  thickness,  the  numerical  grid  size  must  be  quite 
small.  With  a  small  grid  size  and  the  small  steps  in  time 
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needed  for  stability  of  the  solution,  the  numerical 
calculations  for  Hall's  9.3  J  pulse  only  went  out  to  about 
one  third  of  the  pulse  time.  This  is  even  with  the  maximum 
allotted  amount  of  computer  resources.  Nevertheless,  it  is 
felt  that  the  data  obtained  provides  fairly  accurate 
results  since  the  most  intense  laser  energy  is  provided 
very  early  during  the  laser  pulse  (during  the  "spike  of  the 
pulse") . 

In  summary,  this  study  was  directed  towards  examining 
the  effects  of  convection  heat  transfer,  radiation  heat 
transfer,  and  vapor  condensation  on  the  thermal  coupling  of 
a  pulsed,  CO2  laser  to  a  "cold"  Aluminum  surface.  The 
experimental  background  is  Hall's  9.3  J  pulse  (2),  and  the 


II.  The  Experiment 


The  experimental  background  for  this  research  is  Hall's 
famous  9.3  J  pulse  experiment,  in  which  he  noticed  that  the 
thermal  coupling  of  a  CO2  laser  to  an  Aluminum  target 
increased  from  about  3%,  for  a  cold  target  at  CW 
intensities,  to  about  17.5%,  for  the  same  target  exposed  to 
a  9.3  J  pulse.  The  documentation  for  this  experiment  is 
found  in  reference  (2).  This  reference  lists  the  results 
of  several  experiments,  one  of  which  is  the  thermal 
coupling  measurement  being  studied  here.  The  experiments 
were  performed  at  the  Boeing  Aerospace  Laboratory  under 
contract  to  the  the  United  States  Air  Force  Weapons 
Laboratory.  The  laser  used  for  the  experiments  was  a  Marx 
Double  Bank  CO2  laser  of  10.6  micron  wavelength 
radiation,  operated  in  a  pulsed  mode,  with  a  focal  length 
of  28  centimeters,  and  a  focused  spot  radius  of  0.126 
centimeters  at  the  target  surface.  The  irradiated  target 
was  a  2024-T3  Aluminum  alloy  disc  of  5  centimeters 
diameter,  and  a  thickness  of  0.028  centimeters. 

The  back  surface  temperature  of  the  target  was 
monitored  with  a  0.005  inch  diameter  chromel-alumel 
thermocouple  spot-welded  to  the  target.  From  this  measured 
data,  Hall  was  able  to  find  the  absorbed  energy  density  of 
the  target  as  a  function  of  radial  position  from  the  center 
of  the  irradiated  area  (Fig.  1) .  By  integrating  over  the 
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FIG.  is  HALL'S  ABSORBED  ENERGY  DENSITY  PROFILE 


area,  he  determined  that  the  total  laser  energy  residing  in 
the  target  was  1.63  J,  or  17.5%  of  the  incident  beam 
energy,  9.3  J.  By  following  the  above  procedure,  one 
includes  in  a  thermal  coupling  calculation  only  that 
portion  of  laser  energy  which  resides  in  the  target  after 
irradiation,  such  that  the  temperature  of  the  target  is 
increased.  As  such,  impulsive  forces  delivered  to  the 
target  by  pressure  relaxation  are  not  included  in  a  thermal 
coupling  calculation;  rather,  they  are  referred  to  as 
mechanical  coupling  phenomena. 

In  order  to  duplicate  Hall's  data  numerically,  it  is 
necessary  to  know:  the  temporal  variation  of  the  beam,  the 
spatial  variation  of  the  beam,  and  the  peak  flux,  or 
intensity,  of  the  beam.  Hall  presents  a  plot  of  the 
normalized  power*  of  the  laser  used  versus  time.  The 
hydrocode,  as  explained  later,  needs  to  know  the  normalized 
flux  as  a  function  of  time  but,  since  flux  is  just  power 
per  unit  area,  the  two  functions  are  the  same.  To  obtain  a 
functional  expression  for  the  normalized  flux.  Hall's  plot 
was  traced  onto  graph  paper  and  the  resulting  curve  was  fit 
using  the  method  of  least  squares  as  given  by  Book  (13). 

The  following  curve  fit  was  obtained  and,  though  it  is 
tedious,  it  agrees  almost  exactly  with  Hall's  plot: 


*  Normalized  to  the  peak  value. 


p(t)/p*  = 


0.07339  +  1 . 6579t  -  1.2480t2 

+  0.3845t3, 

0.8483  +  0 . 1765t  -  0.05158t2, 

1.5577  -  0.2599  +  0.01508t2  , 

-0 . 025t  +  0.6432 

0 


0  <  t  <  1.5 

1.5  <  T  <  3 

3  <  t  <  8  (1 

8  <  t  <  25.73 

t  >  25.73 


where  p*  is  the  maximum  power  in  the  pulse,  t  is  in 
microseconds,  and  the  pulse  time  is  25.73  microseconds. 
Equation  (1)  is  plotted  in  Fig.  2.  The  remaining  needed 
parameters,  spatial  beam  variation  and  peak  flux,  are  not 
given  explicitly  by  Hall,  though  he  does  say  that  the  beam 
is  Gaussian  with  a  spot  radius  of  0.126  centimeters,  and 
that  the  total  energy  in  the  beam  is  9.3  J. 

The  peak  power,  p* ,  may  be  found  from  equation  (1)  and 
the  total  beam  energy.  Since  energy  is  power  integrated 
over  time. 


E 


9.3  J 


p(t) / p*dt , 


(2 


where  Tp  is  the  pulse  time  of  the  laser,  given  above.  By 
solving  equation  (2) ,  p*  was  found  to  be  0.976121  MW. 
Also,  the  energy  in  the  beam  at  time  t  is  given  by: 


NORMALIZED  POWER 


E  ( t)  =  p* 


(3) 


J  P(  T  )  /p*d  T  , 
0 


where  T  is  now  only  a  dummy  integration  variable.  Equation 
(3)  is  plotted  in  Fig.  3. 

As  stated  above,  the  beam  is  Gaussian  with  a  spot 
radius  of  0.126  centimeters.  The  spatial  form  of  such  a 
beam  is  given  by  Ready  [12}  to  be: 

ENV(r)  =  exp(-ar2),  (4) 

where  ENV(r)  is  the  ratio  of  the  flux  at  radius,  r,  to  its 
peak  value  at  the  center  of  the  beam  (r  *  0) ,  and  a  is  a 
constant  related  to  the  spot  radius.  The  spot  radius  is 
known,  but  there  are  two  common  definitions  of  the  spot 
radius  for  a  Gaussian  beam:  (a)  the  radius  at  which  the 
flux  is  equal  to  exp(-l)  times  its  peak  value  and;  (b)  the 
radius  at  which  the  flux  is  equal  to  exp (-2)  times  its  peak 
value.  It  is  necessary,  therefore,  to  determine  which 
definition  Hall  uses  in  his  report.  He  states  that  the 
Gaussian  width  of  his  calculated  absorbed  energy  density 
profile  is  4.38  millimeters,  and  by  referring  to  Fig.  1, 
one  sees  that  this  radius  corresponds  to  definition  (a) 
above.  It  is  assumed,  then,  that  the  same  definition 
applies  to  the  beam  spot  radius.  From  equation  (4)  and 
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definition  (a) , 


ENV (0.126)  =  exp  -{a (0 . 126) 2}=  exp(-l). 

Thus,  a  =  1/(0.126)2  =  62.988 »63.  The  spatial  variation 
of  the  beam  becomes , 


ENV (r)  =  exp (-63r2)  .  (5) 

Equation  (5)  is  plotted  in  Fig.  4. 

The  only  beam  parameter,  needed  by  the  hydrocode,  which 
remains  to  be  found  is  the  peak  laser  flux.  As  mentioned 
previously,  this  value  is  not  given  by  Hall.  It  is  found 
by  assuring  that  the  laser  flux,  integrated  over  space  and 
time,  is  equal  to  the  total  beam  energy  of  9.3  J.  At  a 
given  time,  the  flux  incident  on  the  mesh  is  given  by 
equation  (1)  to  be, 


I(t)  =  IMAXp ( t ) /p* ,  (6) 

where  IMAX  is  the  peak  laser  intensity.  From  equation  (5) , 
the  flux  varies  spacially  as: 


I (r,t)  =  IMAX ( p ( t) / p* )  exp { -63r 2 ) . 


(7) 


The  integral  of  equation  (7),  over  all  space  and  time,  must 


equal  the  total  beam  energy  of  9.3  J.  Therefore, 


I(r,t)  d 9  dtrdr  =  9.3  J. 


(8) 


Since  I(r,t)  is  not  a  function  of  the  angle,  $,  and  IMAX 
is  a  constant,  equation  (8)  reduces  to: 


2 It  IMAX 


oo  *^p 

f  S  p(t)/p*  r  exp(-63r2)  dtdr 
0  0 


9.3  J. 


(9) 


By  solving  equation  (9),  IMAX  is  found  to  be  19.571 
megawatts  per  square  centimeter,  which  is  approximately  20 
megawatts  per  square  centimeter. 

In  summary,  the  hydrocode  needs  to  know  three  beam 
parameters  in  order  to  determine  the  laser  energy  deposited 
above  the  target.  It  needs  to  know  the  peak  laser  flux, 
the  temporal  variation  of  the  beam,  and  the  radial 
variation  of  the  beam.  Should  the  reader  be  concerned,  the 
axial  variation  of  the  beam  is  determined  by  the  hydrocode, 
as  explained  in  Section  III.  In  this  section,  the  above 
parameters  were  derived  from  Hall’s  report.  The  temporal 
variation  is  given  by  equation  (1),  the  spatial  variation 
is  given  by  equation  (5) ,  and  the  peak  flux  is 


Ill .  The  Hydrocode 

It  was  stated  in  the  introduction  that  the  numerical 
calculations  needed  for  this  study  are  performed  by  a 
three-dimensional  hydrocode.  Henceforth,  the  hydrocode 
shall  be  known  simply  as  HYDRO.  In  this  section,  HYDRO  is 
explained  in  detail  and  all  of  its  equations  are  derived. 
The  version  of  HYDRO  used  (and  modified)  by  the  author  was 
obtained  from  Jackson*.  It  differs  some  from  the  original 
version  which  is  briefly  explained  by  Nielsen  (111. 

HYDRO  solves  the  unsteady  inviscid  flow  equations, 
yielding  flow  properties  outside  the  boundary  layer.  The 
laser  beam  is  assumed  to  be  cylindrical  and  the  flow  is 
assumed  to  be  axisymmetric  about  the  beam  axis,  as  shown  in 
Fig.  3.1. 


/ 

/ 

Target  / 
/ 
/ 
/ 


Laser  Flux 
I (r ,t) 


Fig  3.1  Cylindrical  laser  beam  impinging  on  a 
metal  target. 


The  equations  solved  are  the  "control  volume"  equations  for 
the  conservation  of  mass,  momentum,  and  energy.  In 
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deriving  these  equations,  the  control  volume  of  interest  is 
a  cylindrical  shell,  as  shown  in  Fig.  3.2. 


(*+■  dr )  d  # 


Fig.  3.2  Control  volume  used  for  derivation  of  the 
hydrodynamic  flow  equations. 


Following  the  method  outlined  by  Reynolds  and  Perkins 
[30},  the  equations  may  be  derived  in  the  following  manner: 


(a)  Mass: 


£mg  -  2  Mi  =  - 


dM/d  t)cv,  where 


the  subscript,  i,  means  "in" 
the  subscipt,  o,  means  "out", 

is  the  net  mass  flow  rate,  and 
a  M/a  t)cv  is  the  time  rate-of-change  of 


mass  in  the  control  volume. 


(b)  Momentum:  Z^M0  -  +  dM/dt)cv  = 


£f, 


where 


the  subscript,  i,  means  "in", 
the  subscript,  o,  means  "out", 

Lh  is  the  net  linear  momentum  flow  rate, 
is  the  net  external  force,  and 
Sh/9  t)cv  is  the  time  rate-of-change  of 
linear  momentum  in  the  control  volume. 


(c)  Energy: 


5Zwi  +  (h+ke+pe)i 

=  2rfW0  +  2q  o  +  Sft0(h+ke+pe)0  + 
dE/d  t )  cv ,  where 


the  subscript,  i,  means  "in", 
the  subscript,  o,  means  "out", 

W  is  the  net  rate  of  work  done  by 
external  forces, 

iL* Q  is  the  net  rate  of  heat  input  (or 
output) , 

h  is  the  mass  specific  enthalpy, 
ke  is  the  mass  specific  kinetic  energy, 
pe  is  the  mass  specific  potential  energy,  and 
dE/d  t)  cv  is  the  time  rate-of-change  of 
energy  in  the  control  volume. 


As  stated  previously,  the  flow  is  assumed  to  be 
axisymmetric  about  the  beam  axis;  therefore,  there  is  no 
angular  (0)  dependence  in  the  hydrodynamic  equations. 

Using  the  manner  outlined  above,  the  flow  equations  are  now 
derived. 


Conversation  of  Mass: 


Em0  "  =  -  dM/3  t)cv. 


(la) 


Referring  to  Fig.  3.2, 


p  urrd  0dz  +  P  u2rd  0dr. 


(2a) 


Using  a  first  order  Taylor  series  expansion  about  the  inlet 
conditions , 


=  p  urrd  0  dz 
p  uzrd  0  dr  + 


+  d/d  r  (p  urrd  0dz)dr 
d/d  z  (  puzrd  0dr)dz. 


+ 


(3a) 


Since  mass  is  density  times  volume, 

dVl/d  t)cv  =  d/d  t  (  p  rdrd  0dz)  .  (4a) 


Substituting  equations  (2a)  through  (4a)  into  equation  (la) 
yields : 


1/r  C  d/dz  (p  rur)  +  d/d  z  (p  ruz)  }  =-  ap/  5t.  (5a) 


Now,  the  axisymetric  divergence  operator,  in  cylindrical 
coordinates  is: 


div  (X)  =  1/r  d/d  r  (rXr)  +  d/d  z  (x2) 


With  this  identity,  equation  (5a)  becomes: 


a  t  ®  -div  ( p  v) 


where  the  velocity,  v,  equals  urer  +  uzez, 


Conservation  of  Momentum: 


Em0  -  +  a  M/a  t)  cv  =  2f. 


Consider  first  the  linear  momentum  in  the  radial  direction. 


Mi  =  53Mur  =  pur2rd  0dz  +  puruzrd  0dz. 


Likewise, 


jL/  M0  =  pur2rd  0dz  +  a/a  r  (p  ur2rd  0dz)dr  + 

p  uruzrd  0dr  +  d/d  z  (p  uruzrd  0dr)dz, 


s 


0M/dt  )cv  =  d/dt(Mur)  =  0/ 0  t  (  p  urrdrd  0dz)  .  (4b) 


Finally,  for  no  body  forces, 


2f  =  prd  0dz  -  (p  +  dp/drdr)  rd  0dz.  (5b) 


Substituting  equations  (2b)  through  (5b)  into  equation  (lb) 
yields : 

d(pu  r)/  d  t  =  -1/r  C  d/d  r(  p  rur2)  +  d/d  z(pruruz)  } 

-  d  p/d  r .  (6b) 


Using  the  definition  of  the  divergence  operator,  equation 
(6b)  becomes: 


d(pur) / d  t  =  -div( p  urv)  -  dp/d  r. 


Similarly,  consideration  of  the  linear  momentum  in  the 
axial  direction  yields: 


0(puz)/flt  =  -div(puzv)  -  5p/5  z. 


Conservation  of  Energy: 


+  (h+ke+pe)  i  =  ]E*fW0  +  iLrfQ0  + 

£mq  (h+ke+pe)  0  +  dE/3  t)cv. 


(lc) 


Since  there  is  no  shaft  work  and  no  viscous  work,  the  work 
terras  are  zero.  Assuming  that  the  change  in  potential 
energy  is  very  small, 

Q  =  dE/d  t)cv  +  X/M0(h+ke)0  -  2£Mi(h+ke)i,  (2c) 

where  Q  is  the  net  heat  input  per  unit  time. 

(h+ke)  i  =  purrd  0dz(h+v2/2)  +  puzrd  0dr  (h+v2/2)  .  (3c) 

Using  a  first  order  Taylor  series  expansion  about  the  inlet 
conditions , 


5jM0(h+ke)0  =  purrd  0dz  (h+v2/2)  + 

d/ d  r  [ p urrd0dz (h+v2/2) ]dr  + p  uzrd  6 dr (h+v2/ 2)  + 
d/ d  z  [ p uzrd  0dr (h+v2/2) ]dz.  (4c) 


Now, 


dE/dt  )cv  =  d/dt  (  p  (e+v2/2)  rdrd  0  dz]  ,  (5c) 


where  e  is  the  mass  specific  internal  energy.  Substituting 
equations  (2c)  through  (5c)  into  equation  (lc)  yields: 


Q/rdrd  0dz  =  1/r  d/d  r  [  p  rur  (h+v2/2)  ]  +  (6c) 

1/ra/dz  [pruz(h+v2/2)]  +d/dt  [ p  (e+v2/2)] . 


Now,  rdrd  0 dz  is  the  volume  and  if  Q  is  the  net  heat  input 
per  unit  mass,  per  unit  time. 


Q/rdrd  0dz  =  p  Q. 


To  be  consistent  with  Nielsen's  [11}  documentation,  the 
following  definition  is  made: 


E  =  p(e+v2/2) . 


By  definition,  h  =  e  +  p//>  ;  therefore, 


h  +  v2/2  =  1/p  (E+p) 


Substituting  equations  (7c)  through  (9c)  into  equation  (6c) 
yields : 


d  E/at  =  -l/r  d/dr [rur (E+p) ]  -  d/dz [uz (E+p) ]  +  p  Q.  (10c) 


Using  the  definition  of  the  divergence  operator,  equation 
(10c)  becomes: 


aE/dt  =  -div  [(E  +  p)  v]  +PQ. 


(13) 


The  "control  volume"  conservation  equations  are: 


dp/dt  =  -div  (  p  v)  , 

5(Pur)/3t  =  -div  (ur  p  v)  -  dp/  dr, 
d( P  uz) /  d  t  =  -div (uz  p v)  -  dp/ d  z, 
dE/d  t  =  -div  [(E+p)v]  +  P  Q. 

There  are  only  four  equations  in  the  five  unknowns,  p  , 
ur,  uz,  E,  and  p.  An  additional  equation  is  needed  in 
order  to  simultaneously  find  the  five  unknowns.  The 
equation  used  is  the  perfect  gas  equation  of  state  for  air, 
p  =  P  RT.  Assuming  that  the  gas  is  calorically  perfect,  R 
=  Cp  -  cv  =  cv( y  -1) ,  where  Y  is  the  ratio  of 
specific  heats.  Also,  for  the  same  assumptions,  e  *  cvT; 
therefore,  p  =  Pe(7-1).  Since  E  =  P(e+v2/2) , 

p  =  (E  -  v2/2)  (  y-1)  ,  (14) 

where  v  =  (ur2  +  uz2)^2.  Equation  (14)  augments 
equations  (10)  through  (13).  The  hydrodynamic  equations 
solved  by  HYDRO  are  thus  derived. 

Equations  (10)  through  (14)  are  solved  by  HYDRO,  using 
finite  difference  techniques.  Before  the  difference 
equations  are  derived,  a  suitable  mesh  or  grid  must  be 
established.  The  mesh  used  by  HYDRO  corresponds  with  the 
control  volume  shown  in  Fig.  3.2.  The  axial  mesh  is  set  up 


(10) 

(11) 

(12) 

(13) 


such  that  the  axial  cells  have  a  constant  thickness,  Az, 
in  the  axial  direction,  as  shown  in  Fig.  3.3  below.  The 
radial  mesh  is  set  up  such  that  the  cylindrical  cells  have 
a  constant  thickness,  At  ,  in  the  radial  direction,  as 
shown  in  Fig.  3.4  below.  For  simplicity  Az  -  At  =  Ax. 

zl.  2  2 


Target 

Ax/  2  _ 

,  A  x  > 

Zi  =  (i-1/2)  Ax 

Fig.  3.3  Schematic  of  axial  mesh. 


rj  =  (j-1)  Ax, 

r^  corresponds  to  the 
center  of  the  beam. 


Fig.  3.4  Schematic  of  radial  mesh. 


To  derive  the  difference  equations  for  a  flux-divergent 
quantity  q,  such  that  dq/dt  =  -divf.  Gauss'  divergence 


theorem  is  used.  The  theorem  is  stated  as: 


JJJ 


divf  dV  = 


For  a  differential  volume  element,  the  theorem  reduces  to: 


divf  (V)  =  f  0A0  -  fj.Ai  where 

V  is  the  volume  of  the  element, 

f0A0  is  the  amount  of  q  leaving  the 
element  per  unit  time,  and 

f^A^  is  the  amount  of  q  entering  the 
element  per  unit  time. 

Consider  first  the  divergence  of  a  quantity,  q,  in  the 
radial  direction.  The  volume  element  is  a  cylindrical 
shell  of  height  Ax,  inner  radius  r,  and  outer  radius 
r  +  Ax.  By  Gauss'  theorem, 

divf  =  f  (r+  4  x)  2?r(r+  Ax)  Ax  -  f  (r)  2n  r  A  x 
[  7T(r+  Ax)2  -  nr2]  Ax 

This  equation  reduces  to: 

divf  =  1  /Ax  (1+  A  x/2r)“l[(l+  Ax/ r)  f(r+  4x)-f(r)]. 


For  a  flux-divergent  quantity,  q,  dq/dt  =  -div£;  therefore 


dq/dt  =  -1/4  x  [1  +  4  x/(2r))*1[(l+4  x/r)f  (r+4x)-f  (r)  ]  .  (15) 


Consider  now  the  divergence  of  q  in  the  axial  direction. 
This  case  is  simpler  since  V  =  a4x,  and  Gauss'  theorem 
yields ; 


dq/dt  =  -1/4  x  [f  (z+ 4  x) -f  (z)  ]  .  (16) 

Equations  (15)  and  (16)  may  be  combined  into  a  single, 
equation  by  defining  a  parameter  1,  such  that  1=1 
corresponds  to  the  axial  direction  and  1=2  corresponds  to 
the  radial  direction.  The  combined  equation  is: 

dqi/dt  =  -1/4  x[l+((l-l)  4x)/ (2xx)  I”3-  X 

Ul+Ul-1)  4  x) /xi  ]f  <xi  +  4x)-f  (xx) }.  (17) 

Equation  (17)  does  not,  however,  apply  to  the  radial 
flux  of  quantities  from  the  first  radial  cell.  The  center 
of  this  cell  is  at  the  center  of  the  beam?  therefore,  there 
is  only  one  boundary  through  which  quantities  are 
transported.  For  this  cell.  Gauss'  theorem  yields; 

divf  7r(  4x)2/4  =  JT4xf3/2-  (18) 

Equations  (17)  and  (18)  are  used  by  the  difference  method 
to  update  hydrodynamic  quantities  as  they  vary  with  both 
time  and  spatial  position  in  the  mesh. 


The  numerical  method  used  to  update  quantities  is  the 
Lax-Wendroff  two-step  technique,  which  is  accurate  to  order 
(  At/  4x)2.  A  detailed  explanation  of  this  method  may  be 
found  in  reference  {14}.  Here,  the  method  is  explained 
only  to  the  extent  that  the  interested  reader  will  be  able 
to  understand  how  HYDRO  does  the  differencing.  The  two 
steps  employed  are: 

(1)  qt  +  1/2  =  qt  +  dq  /dt(dt/2),  and 

n  +  1/2  n  +  1/2  n 

(2)  qt+1  =  qt  +  dqn-i/2/dt  {  A  t)  . 

n  n 

In  these  equations,  q*  is  the  value  of  a  flux-divergent 

n 

quantity  at  spatial  position  n  and  time  step  (cycle)  t,  and 
ft  is  the  corresponding  flux  component  in  the 
direction  of  interest.  Also,  qn+i/2  defined  to  be 
the  arithmetic  mean  of  qt  and  qn+i*  Fig.  3.5  on 
the  following  page  will  help  the  reader  visualize  the 
notation  used  above. 


Fig.  3.5  Schematic  of  Lax-Wendroff  two-step  technique 


Without  going  through  the  algebra,  the  two  Lax-Wendroff 
steps,  with  dqn/dt  given  by  equations  (17)  and  (18), 
become : 


qt-1/2  =  (qt  +  qt ) / 2  - (  A  t/ 2  A  X)  [ 1 +  ( 1-1 ) /  (  2  ( n- 1 )  )  ] " 1 
n+1/2  n  n+1 


C [1+ (1-1) / (n-1) ]ft  -  ft3, 

n+1  n 


( 


qt  +  i=  qfc  -  [  4t/  2ix)  [1+  (1-1)  /  (2n-3)  ]  -!(  [1+  (1-1)  /  (n-3/2)  ] 
n  n 


ft+1/2  _  ft+l/2j 
n+1/2  n-1/ 2 


for  all  cells  except  the  first  radial  cell,  and 


qt^l/2  =  l/2[qt  +  qt  _  (2  At/  A  X) f ( A  x ) 1 ,  ( 

qt+1  =  qt  _  (  4t/  4x)  4ft','1/2f  ( 


for  the  first  radial  cell. 


Equations  (19)  through  (22)  are  used  to  update  density, 
energy  density,  and  the  axial  component  of  the  momentum. 

The  radial  component  of  the  momentum  is  troublesome, 
though,  because  it  contains  a  gradient  term,  -  dp/  dr,  and 
must  be  handled  differently.  This  term  is  handled  by 
assuming  that  dp/d  r  is  approximately  equal  to  4p/  At. 
With  this  assumption,  the  two  Lax-Wendroff  steps  become, 
for  q  =  p ur: 

qt+1/2  =  (qt+l/2)f_d  -4t/(2  4x)  (pfc  -  pfc)  ,  (19a) 

n+1/2  n+1/2  n+1  n 


qt  +  1  =  (qt+l)f  -  At/ A  x(pt+3-/2  -  pt+1/2),  (20a) 

n  n  n+1/2  n-1/2 

=  f-d  -  At/(24x)  (pj  -  pt)  ,  (21a) 

3/2  3/2  21 

qt+1  =  (qt+1)  £_d  -  4t/4x(pt+l/2  -  pt+1/2),  (22a 

1  1  1  3/2  1/2 


In  equations  (19a)  through  (22a),  the  subscript,  "f-d", 
refers  to  the  flux-divergent  term  given  by  equations  (19) 
through  (22)  . 

According  to  Nielsen  (11),  it  is  sometimes  convenient  to 
difference  pressure,  instead  of  energy,  as  an  independent 
variable.  This  is  done  to  avoid  negative  pressures  and 
temperatures  near  the  target  when  laser-induced  boiling  is 
taking  place.  In  pressure-  differencing,  equation  (13)  is 


transformed  by  using  the  perfect  gas  equation  of  state 
(equation  (14)),  the  mass  continuity  equation  (equation 
(10)),  and  Euler's  equation  in  vector  form.  Euler's 
equation  is  stated  as: 


dv/  d  t  =  -(v  •  7)  v  -  7p/p  . 


With  this  transformation,  equation  (13)  is  replaced  by: 


dp/d  t  =  -div(pv)  -  (y-l)pdivv  +  (  y  -1)  pQ. 


To  make  equation  (23)  consistent  with  the  flux-divergent 
form,  dq/dt  =  -divf,  the  following  adaption  is  made: 


ap/dt  =  -div  (P+(  y-i)pc)v  +  (y-i)pQ, 


where  pc  is  a  constant  pressure  evaluated  at  the  center  of 
a  cell,  taken  to  be  the  arithmetic  mean  of  the  pressures  at 
the  boundaries  of  the  cell.  Fig.  3.6  helps  clarify  the 
meaning  of  pc .  One  of  the  advantages  of  the  flux-divergent 
equations  is  that  the  flux  of  a  quantity  out  of  one  cell  is 
the  flux  of  that  quantity  into  the  next  cell.  Due  to  the 
presence  of  pc  in  equation  (24)  ,  one  sees  that  the 
"pressure  flux"  out  of  a  cell  will  not  be  the  "pressure  flux 
into  the  next  cell.  HYDRO  handles  this  problem  by 
multiplying  the  flux  out  of  a  cell  by  a  "correction  factor," 
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such  that  it  approximates  the  flux  into  the  next  cell.  The 
correction  factor  is  described  below. 


cel  ln 

celln+i 

Pc 

Pc' 

*  - » 

• 

Plef  t 

Pright 

vlef  t 

vright 

cell  width  =  4x 


Fig.  3.6  Schematic  of  cell-centered  pressure-differencing. 

If  the  pressure  flux  at  the  left-hand  boundary  and  that  at 
the  right-hand  boundary  are  given  respectively  by, 

fl  =  (Pl  +(  Y-l)pc)vi  and  fr  =  (pr+  (  y-l)Pc)vr, 

one  sees  that  the  flux-divergent  part  of  dp/dts: 

-f  (fr-f  i)  /  Ax.  As  stated  above,  frn  * f ln+1 ' 
therefore,  frn  is  multiplied  by  a  "correction  factor" 
such  that  it  will  approximate  fin+i*  there  are  two  of 
these  correction  factors:  one  for  space-interval  steps  and 
one  for  time-interval  steps.  Referring  to  Fig.  3.6,  it  is 
evident  that  the  space-interval  factor  is  given  by. 


CF  =  (pr+{  y-l)p’c)/(Pr+(  y-l)Pc>- 


(25) 


In  equation  (25),  Pc  =  (pr+Pi)/2  and  pc '  =  (pn+2+Pr)/2» 
where  the  subscripts  r  and  1  mean  "right"  and  "left", 
respectively.  Also  pr  =  pn+i  and  =  pn,  where  n 
is  the  cell  number.  Correction  factors  for  both  the  radial 
and  axial  directions  are  found  by  equation  (25);  however, 
the  time-interval  correction  factor  is  handled  somewhat 
differently.  It  is  defined  by  equation  (25),  but  pc '  is 
now  zero  since  no  "time-pressures"  have  been  found  past 
time  t, 

Pr  =  pt+1/2,  extrapolated  from  pt+1/2,  and 
n+1/2  n 

Pi  =  pt+l/2,  known  from  previous  calculations.  The  time- 
1  n-1/2 

interval  correction  factor  is  thus  given  by, 

CF  =  (pt+1/2)/  pt+l/2  +  0.5 (  y  -1) (pt+l/2  +  pt+1/2).  (26) 

n+1/2  n+1/2  n+1/2  n-1/2 

For  pressure-differencing,  equation  (24)  ,  augmented  with  the 
correction  factors  given  by  equations  (25)  and  (26) , 
replaces  equation  (13)  of  the  hydrodynamic  flow  equations. 

In  order  to  difference  either  energy  or  pressure  with  the 
same  code,  the  variable  NDIF  is  introduced.  If  NDIF  =  0, 
pressure  is  differenced,  and  if  NDIF  =  1,  energy  is 
differenced.  Formulas  which  differ  in  the  two  cases  are 
written  in  the  form, 

X  =  (energy  formula)  x  NDIF  +  (pressure  formula)  x  (1-NDIF) . 


As  stated  previously,  when  a  laser  of  sufficient 
intensity  strikes  a  target  surface,  boiling  and 
vaporization  of  the  target  may  occur.  Consider  a  control 
volume  (as  shown  in  Fig.  3.2)  near  the  target  surface, 
consisting  of  air  and  target  vapor.  If  one  assumes  that 
the  mixture  is  dilute,  the  species  mass  conservation 
equation  for  target  vapor  yields*, 

dpv/3t  =  -div ( p  vvv) .  (27) 

The  development  of  equation  (27)  may  be  found  in  Vincent i  and 
Kruger  (15).  According  to  the  approach  taken  by  Anisimov  C 16  3 
conservation  of  energy  in  the  control  volume  yields, 

vd  =  fb/{ p 0 (L+CTV) ] ,  (28) 

where  vd  is  the  velocity  of  the  vaporization  front  into 
the  target,  fb  is  the  incident  laser  flux,  and  pQ,  L, 

C,  and  Tv  are  the  target  density,  latent  heat  of 
vaporization,  specific  heat,  and  vaporization  temperature. 

To  determine  Tv  self-consistently  at  the  elevated  recoil 
pressure  of  the  blow-off  vapor,  equation  (28)  is  augmented 
by, 

Vd  =  Cs  exp (-L/CTV) ,  (29) 

*  In  this  development,  the  subscript,  v,  denotes  target  vapor. 


where  Cs  is  the  speed  of  sound  in  the  target  material. 
Steady-state  conservation  of  mass  at  the  target  surface 
requires  that  pvvv  =  P ovd'  where  is  given  by 


equation  (28).  Thus  pv  is  found  to  be, 

Pv  =  ffa/  [Vy (L+CTy)  ]  . 

Substitution  of  equation  (29)  into  equation  (28)  yields 

fb /[  Po (L+CTV) ]  =  Cs  exp (-L/ CTV) . 

From  equation  (31) ,  CTV  is  expressed  by  the 
transcendental  equation, 

CTV  =  L/ln  [  p0Cs(L+CTv)/fb] . 

Assuming  that  the  blow-off  vapor  is  a  perfect  gas,  both 
calorically  and  thermally. 


pv  =  (  y— 1)  Pvev  =  PVCTV(  y-1)  • 


Now, 


Pv  =  Pvvv^' 


so  that 


vv  =  (Pv/Pv)l/2  =  (  y-1)  1/2  (CTV)  1/2. 


Substituting  this  result  into  equation  (30)  yields. 


Pv  =  fb/Ifi  (  y-D  1/2  <CTv)1/2], 


whereas  L+CTV.  The  total  energy  in  the  vapor  is, 


Ev  =  Pv(ev+vv 2/2)  =  0.5  PVCTV(  y +1)  .  (3 


Substitution  of  equation  (33)  into  equation  (34)  yields. 


Ev  =  fb(  r-l)/(2  Q  )  [CTv/(  y-l)]l/2.  (35 


The  mass  flux  at  the  surface  is. 


Ms  =  Pvvv  =  fb/ O 


One  sees,  therefore,  that  the  evolution  of  target  vapor  is 
governed  by  the  following  equations: 


Pv  -  fb/tfl  <r-DL/2  (CTv)1/2], 


Ev  =  fb(  y  +D/  <2  Q  )  [CTV/  (  y-1)  ]  1/2, 


With  the  velocity  of  the  vaporization  front  into  the 
target  given  by  equation  (28)  ,  it  is  possible  to  determine 
the  depth  to  which  the  laser  has  penetrated  the  target  at 


time,  t.  Denoting  this  depth  by  CHEW(r,t), 


CHEW(r,t) 


t 

J  vddT 

0 


tfb  (r  ,t)  /  (pQ  Q  )  4  t. 


(36a) 


The  extent  to  which  CHEW  is  advancing  in  time  is  an 
indication  of  how  plasma  effects  are  shielding  the  surface 
from  laser  radiation.  One  can  therefore  determine  the  time 
at  which  a  surface-shielding  LSD  is  formed. 

Two  key  assumptions  are  implicit  in  the  derivation  of 
equations  (33)  through  (36) .  First,  it  is  assumed  that  the 
fraction  of  incident  laser  radiation  (  a  )  absorbed  by  the 
target  is  unity,  although  in  general  the  derivation  of  the 
equations  need  not  have  a  value  of  unity.  According  to 
Schwirzke  C 2 3 } ,  a  critical  density  sublayer  of  plasma  is 
formed  at  the  target  surface  in  the  nanosecond  range  and  the 
mechanism  of  unipolar  arcing  assures  that  a=  1  once  this 
sublayer  is  formed.  The  time  frame  of  the  laser  pulse  under 
study  is  in  the  microsecond  range?  therefore,  a  should 
approach  1  rapidly.  When  a  LSD  is  ignited,  a  actually 
decreases  due  to  dissipation  of  the  critical  density 
sublayer.  However,  the  amount  of  laser  attenuation,  as 
discussed  below,  assures  that  the  flux  incident  on  the 
surface  rapidly  decreases  as  the  LSD  forms.  For  these 


reasons,  it  is  felt  that  the  first  assumption  is  valid  for 
the  problem  under  study.  The  second  assumption  is  that  the 
radial  and  axial  conduction  of  heat  into  the  target  surface 
is  negligible  over  the  time  of  the  pulse.  From  dimensional 
analysis  of  the  one-dimensional  heat  conduction  equation, 
the  diffusion  time  required  for  radial  conduction  may  be 
estimated  as  tr  sr s2/*  ,  where  rs  is  the  laser  spot 
radius  and  K  is  the  thermal  diffusivity  of  the  target 
material.  Hall  £23  gives  rs  and  K  to  be  0.126  cm  and 
0.54  cm2/sec,  respectively.  Thus,  tr  »  30,000 
microseconds  which  is  about  three  orders  of  magnitude 
larger  than  the  laser  pulse  time  of  about  26  microseconds. 
This  may  also  be  compared  to  the  characteristic  time  in  the 
direction  into  the  target.  By  the  same  reasoning  this 
relaxation  time,  tz  s  d2/*  ,  where  d  is  the  target 
thickness  (.028  cm),  is  tz  st  1450  microseconds.  It  is 
felt,  then,  that  the  second  assumption  is  valid. 

As  a  laser  beam  propagates  through  a  gaseous  medium, 
some  of  the  laser  energy  is  deposited  in  the  medium.  The 
amount  of  energy  deposited  depends  on  the  absorption 
coefficient  of  the  medium.  Two  absorption  coefficients  are 
used  by  HYDRO;  a  "cold"  absorption  coefficient  representing 
absorption  in  molecular  bands  of  the  target  vapor,  and  a 
"free-free"  coefficient  representing  inverse  bremsstrahlung 
absorption  by  electrons.  These  coefficients  are  defined  as 


^cold  =  Ci  Pv» 


(37) 


kff  =  C2ne2/(kT)3/2, 

where  ne  is  the  electron  equilibrium  number  density.  The 
constant,  Cj_,  is  chosen  such  that  p  0  is  the 
absorption  coefficient  of  the  bulk  solid.  The  constant, 

C2,  is  given  by  Zeldovich  and  Raizer  (17]  to  be, 

C2  =  4/3(2  7r/3m)1/2e6h2/[mc(h  v  )2]  ss  4.05  x  10"53,  (38) 

where  m  is  the  electronic  mass,  c  is  the  speed  of  light,  h 
is  Planck's  constant,  e  is  the  electronic  charge,  and  vis 
the  laser  frequency.  The  temperature  in  a  cell,  needed  by 
equation  (37)  is  found  using  the  ideal  gas  equation  of 
state,  given  by  p^  =  p^kT/m±  for  species  i,  where 
m^  is  the  molecular  weight  of  the  ith  species.  By 
Daltons  Law  of  partial  Pressures,  p  =  Pi+p2,  where  the 
subscript  1  corresponds  to  air  and  the  subscript  2 
corresponds  to  target  vapor.  Therefore, 

kT  =  mim2p/  (  P^m2+  p2mi)  .  (39) 

Equation  (39)  is  the  equation  of  state  used  by  the  original 
version  of  HYDRO.  In  order  to  account  for  high  temperatures 
and  pressures,  this  equation  was  modified  by  the  author  to 


include  the  compressibility  factor  for  air.  The 
compressibility  factor,  Z,  is  defined  by  p  =  pzkT/m.  With 
the  inclusion  of  Z,  equation  (39)  becomes, 

kT  =  m^m2P/(Pim2Z  +  P2ml)  •  (40) 

Now,  p.  pj_,  and  are  found  from  the  hydrodynamic 

flow  equations.  The  temperature  may  be  found  by  an 
iterative  process,  such  that  equation  (40)  is  satisfied. 
Once  the  temperature  is  found,  the  absorption  coefficients 
given  by  equation  (37)  are  determined. 

Laser  attenuation  and  energy  deposition  are  illustrated 
by  Fig.  3.7  below. 
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Fig.  3.7  Schematic  of  laser  attenuation  and  energy 
deposition . 

th 

In  Fig.  3.7,  fj.  +  i  is  the  laser  flux  entering  the  i 
cell  and  fj.  is  the  flux  leaving  the  cell.  Also,  K  is  the 
total  absorption  coefficient  in  the  cell  and  E  is  the  rate 
of  energy  deposition  in  the  cell.  Using  a  first  order 


Taylor  series  expansion  about  the  cell  entrance,  the  flux 
leaving  an  arbitrary  cell  may  be  expressed  by. 


f ou t  =  f  +  (  df/dx)dx,  (41) 

where  f  is  the  flux  entering  the  cell.  Now,  part  of  the 
flux,  f,  passes  through  the  cell  unattenuated;  whereas,  part 
of  it  is  absorbed  while  in  the  cell.  If  the  amount  absorbed 
is  given  by  K^fdx,  then  the  following  equation  must  be 
true : 


fout  +  Kifdx=  f.  (42) 

Substituting  equation  (41)  into  equation  (42),  one  sees  that 

f  =  f  +  (  6  f/dx)dx  +  K^fdx.  (43) 

The  f's  in  equation  (43)  cancel  each  other  and  the  resulting 
differential  equation  is  separable.  By  separating  the 
variables,  one  obtains, 


df/f  =  -K^dx  (44) 

Equation  (44)  may  be  integrated  from  the  inlet  condition, 
cell  i+1,  to  the  exit  condition,  cell  i.  The  result  of  the 
integration,  when  rearranged  is. 


f  i  =  f  i  +  iexp  (-Ki  Ax)  .  (45) 

Thus,  the  amount  of  laser  attenuation  is  determined  by 
equation  (45)  ,  where  is  taken  to  be  the  sum  of  Keff 
and  Kcoid/  given  by  equation  (37)  .  With  f^  known,  the 
time  rate-of-change  of  energy  density  in  the  cell,  due  to 
laser  deposition,  is  simply  df/  d  x,  which  is  approximately 
(f i+l-f i) /  A  x.  Therefore, 

Ei  =  (f i+1-f i) /  A  x.  (46) 

In  order  to  solve  equations  (45)  and  (46)  ,  the 
equilibrium  electron  number  density,  ne,  must  be 
calculated.  The  calculation  of  ne  is  done  using  a  two 
material  Saha  equation  that  allows  only  single  ionization. 
For  single  ionization  only,  there  are  five  species  present 
in  equilibrium:  A,  B,  A+ ,  B+ ,  and  e”.  The  corresponding 
concentrations  are  NA,  Ng,  NA,  Ng,  AND  ne.  Note 
that,  as  the  final  result  of  this  development  will  show,  it 
does  not  matter  which  species  (air  or  target  vapor)  is 
denoted  by  A  or  B.  The  local  equilibrium  equations  are: 

NAneRi(A)  =NAne2RR(A),  (47) 

NBneRI(B)  =Ngne2RR(B),  (48) 


where  Rj  is  he  ionization  rate  and  RR  is  the 
recombination  rate.  Equations  (47)  and  (48)  are  subject  to 
the  "atom  conservation”  equations,  given  by: 


NA  +  NA  =  nA, 


(49) 


Nb  +  nb  =  ng. 


(50) 


NA  +  NA  =  n 


e ' 


(51) 


where  n^  is  the  number  density  of  the  i^  species. 
Substitution  of  equations  (49)  and  (50)  into  equations  (47) 
and  (48)  yields: 

nA  -  NA  -  NAneRR(A) /Rj (A) ,  (52) 

nB  -  Ng  =  NBneRR(B)/Rj (B) ,  (53) 


Addition  of  equations  (52)  and  (53)  yields: 

nA  +  nB  -  (Na  +  NB)  =  ne(NARR(A) /RX (A)  +  NBRR(B) /Rj (B) ] .  (5 

By  substituting  equation  (51)  into  equation  (54),  one 
obtains , 

n  -  ne  =  ne[NARR(A) /RX (A)  +  NBRR (B) /Rx (B) ]  (55) 
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where  n  =  nA  +  ng.  Equation  (55),  expressed  in  terms 
of  nA  and  ng  becomes : 


n  -  ne  =  ne [ [nA(ne+Sg)  +  nB (ne+SA) ] / (ne+SA) (ne+Sg) ] ,  ( 


where  Si  is  defined  as  Rj/Rr  for  the  ith  species. 
Equation  (56)  is  a  cubic  equation  in  ne.  By  letting 
ne  =  x  -  (SA  +  Sg)/3,  one  arrives  at  the  following 
equation  for  x: 


x3  +  ax  +  b  =  0 


a  =  SASB  “  nAsA  "  nBsB  ~  i/3(SA  +  Sg)2, 


b  =  1/27  [2(Sa  +  Sg) 3  -  9(Sa  +  Sg) (SASg  -  nASA  -  ngSg) ] 
-  nSASg. 


Before  equations  (57)  through  (59)  can  be  solved,  SA  and 
Sg  must  be  found.  According  to  Zeldovich  and  Raizer 
(173,  they  are  determined  by  the  following  relations: 


SA  =  G(kT)3/2  exp (-1 A/kT) , 


SB  =  G(kT)3/2  exp (- Ig/kT) , 


( 


where  1^  and  I3  are  the  ionization  potentials  of 
species  A  and  B,  and  the  constant,  G,  has  the  value,  2.98  x 
1()39,  when  kT  is  expressed  in  ergs.  With  and  S3 
determined  by  equations  (60)  and  (61)  ,  the  solution  of 
equation  (57)  has  been  found  by  trial  and  error  to  be: 

x  =  2(-a/3)1/2  cos(0/3),  (61) 

where 

cos  -b/2a(27/-a) i/2.  (62) 

With  x  given  by  equations  (61)  and  (62) ,  ne  is  determined 
by  the  relation,  ne  =  x  -  (S^  +  Sb)/3.  The  free-free 
absorption  coefficient  is  determined  by  equation  (37) ,  and 
once  Kff  is  known,  the  energy  deposition  rate  is  given  by 
equation  (46)  .  The  total  amount  of  laser  energy  deposition 
in  the  mesh  is  therefore  approximated  by: 

E  =  £  E(fi  +  i  -  f i )  4t/4x(Vi),  (63) 

i  t 

where  the  first  summation  is  over  all  axial  cells,  the 
second  summation  is  over  time,  and  is  the  volume  of  the 
ith  cell. 

HYDRO  is  equipped  to  compute  the  mechanical  effects  of  a 
LSD  on  the  surface,  and  even  though  they  are  not  of  interest 


in  this  study,  they  are  discussed  briefly  in  order  that  the 
reader  may  be  able  to  understand  how  HYDRO  operates. 

First,  the  total  impulse  delivered  to  the  target  may  be 
approximated  by: 

I(t)x:  2  (Pj  ~  ^ambient)  Aj  4  t ,  (64) 

t  j 

where  pj  is  the  pressure  in  the  jth  radial  cell  at  the 
target  surface  and  Aj  is  the  corresponding  area  of  the 
cell.  In  equation  (64),  the  first  summation  is  over  time 
and  the  second  summation  is  over  all  radial  cells  at  the 
surface.  The  spot  impulse  is  defined  as  the  impulse 
delivered  over  the  irradiated  area  only  and  is  found  by 
equation  (64)  ,  except  that  the  second  summation  is  now  over 
all  radial  cells  out  to  the  cell  at  which  laser  radiation 
is  effectively  zero. 

The  time  increment.  At,  has  been  mentioned  several 
times  thus  far,  but  it  has  not  been  defined.  It  is 
important  to  note  that,  due  to  the  presence  of  a  LSD,  At 
cannot  be  chosen  arbitrarily  or  the  LAX-Wendroff  technique 
becomes  unstable;  rather,  it  must  be  chosen  such  that 
stability  of  the  solution  is  assured.  This  is  done  by 
locating  the  maximum  mass  specific  energy  in  the  mesh  and 
defining  At  by. 


At  = 


0  Ax/  (E/  p  )1^2, 


(65) 


where  E/p  is  defined  by  equation  (8c)  .  The  factor,  fi  , 
known  as  the  "courant  number"  is  normally  taken  to  be  0.1. 
One  sees,  therefore,  that  equation  (65)  assures  that  a 
disturbance  propagates  no  further  than  one  tenth  of  a  cell 
size  during  At.  For  further  discussion  of  accuracy  and 
stability  criteria,  the  reader  is  referred  to  reference 
(14). 

There  are,  however,  two  types  of  instabilities  which 
will  be  discussed  here.  The  first  arises  when  the  pressure 
gradient  and  velocity  are  in  the  same  direction.  This 
situation  may  occur  near  the  surface  in  the  presence  of 
highly  ionized  target  vapor.  HYDRO  is  modified  to  handle 
this  problem  and  the  modification  is  discussed  in  detail  by 
Nielsen  (11).  It  suffices  to  say  here  that  the  effect  of 
the  modification  is  to  prevent  the  Lax-Wendroff  technique 
from  giving  negative  pressures  and  temperatures  in  the  cell 
previous  to  the  LSD  front.  The  other  instability  arises 
due  to  the  fact  that  HYDRO  is  inviscid  in  nature,  and  no 
dissipative  mechanisms  are  used  in  solving  the  flow 
equations.  With  this  approach,  the  LSD  front  is  a 
numerical  discontinuity,  and  modification  in  the 
neighborhood  of  the  front  is  necessary.  HYDRO  does  the 
modification  by  introducing  an  "artificial  diffusion" 
factor,  discussed  in  detail  in  reference  (11).  Briefly, 
stated,  the  diffusion  factor  has  the  effect  of  spreading 
the  LSD  front  over  several  cells,  while  assuring  that  the 


Rankine-Hugoniot  equations  £18]  are  satisfied  on  either 
side  of  these  cells. 

In  summary,  HYDRO  solves  the  inviscid  flow  equations 
for  a  mixture  of  air  and  target  vapor,  through  which  laser 
energy  is  propagating.  In  this  section,  these  equations 
were  derived  for  the  cylindrical  control  volume  of 
interest.  The  finite-difference  technique  used  to  solve 
the  equations  was  discussed  and  the  difference  equations 
were  derived.  Finally,  the  topics  of  laser-target 
interaction  and  laser  attenuation  were  discussed. 


Convection  Heat  Transfer 


As  stated  in  the  introduction,  the  first  object  of  this 
study  was  to  perforin  detailed  calculations  of  the 
convection  heat  transfer  taking  place  at  the  target 
surface,  under  the  presence  of  a  LSD.  In  this  section,  the 
method  used  to  perform  these  calculations  is  discussed  in 
detail.  Numerous  relations  exist,  by  which  one  can  compute 
the  amount  of  convective  heat  transfer  to  a  surface,  over 
which  a  fluid  is  flowing.  The  relation  that  one  chooses  to 
use  depends  on  the  particular  flow  situation.  A  flow  field 
may  be  considered  to  consist  of  an  inviscid,  far  region,  or 
external  flow  field,  and  a  viscous,  near  region,  or 
boundary  layer.  Resistance  to  the  transfer  of  heat  between 
a  surface  and  a  fluid  is  confined  to  the  boundary  layer. 

In  this  region  the  physical  processes  are  governed  by  the 
viscous  flow  and  energy  equations. 

The  external  flow  field  is  provided  by  HYDRO,  as 
described  in  Section  III.  For  the  problem  at  hand,  the 
field  closely  resembles  that  of  axisymmetric 
three-dimensional  stagnation  flow.  As  given  by  Schlichting 
119),  the  radial  velocity  component  for  such  a  flow  is. 


ur  =  ar. 


where  a  is  a  constant  and  r  is  the  radius  from  the  axis  of 


symmetry.  The  radial  velocity  is  plotted  in  Fig.  5,  from 
which  one  sees  that  the  flow,  out  to  the  LSD  front,  is 
closely  approximated  by  equation  (66)  .  There  exists  an 
exact  solution  to  the  Navier-Stokes  equations  for  3-D 
stagnation  flow.  This  solution  is  given  by  Schlichting, 
but  is  rederived  here  for  completeness.  The  Navier-Stokes 
equations  for  axisymmetric  flow  are: 

u  du /  dr  +  wdu/dz  =  (-1  /p  )  9p/3  r  + 

v[  d2u/  dr2  +  (1/r)  du/  dr  -  u/r2  +  d2u/  dz2] ,  (67) 

u  dw/  d  r  +  wd  w/  d  z  -  (-1/p  )  dp/ d  z  ♦ 

v[d2 w/  3r2  +  (1/r)  dw/dr  +  d2 w/  dz2] ,  (68) 

du/dr  +  u/r  +  dw/d  z  =  0  (69) 

where  v  +  uej-  +  wez.  The  boundary  conditions  are: 

at  z  =  0,  u  =  0  and  w  =  0,  (67a) 

at  z  =  oo,  u  =  U,  (68a) 

where  U  is  the  velocity  of  the  external  flow  field.  Note 
that  equations  (67)  through  (69)  are  for  steady-state, 
incompressible  flow  with  no  body  forces.  Neglecting  the 
effect  of  gravity,  there  are  no  body  forces  for  the  problem 
under  study.  The  assumptions  of  steady,  incompressible 
flow  are,  strictly  speaking,  incorrect  for  this  problem; 
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however,  it  is  felt  that  the  information  provided  with 
these  assumptions  will  be  qualitatively  correct  and  quite 
useful  in  the  analysis  which  follows.  These  ideas  are 
discussed  in  more  detail  later. 

In  ideal  3-D  stagnation  flow, the  stagnation  pressure  is 
related  to  the  local  pressure  by, 

Po  =  P  +  Pi 2(v2)  =  p  +  p/2(a2)(r2  +  4z2)  ,  (70) 

where  pQ  is  the  pressure  at  the  stagnation  point,  p  is 
the  local  pressure,  z  is  the  axial  distance  from  the 
stagnation  point  opposite  the  direction  of  the  impinging 
flow,  and  the  external  velocity,  U  =  arer  -  2azez. 

Here,  a  and  r  are  the  same  as  in  equation  (66).  For  the 
viscous  boundary  layer  flow,  the  following  forms  of  the 
solutions  for  the  velocity  and  pressure  distributions  are 
assumed: 


u 

w 


Po  "  P 


rf *  (z)  , 

(71) 

-2f (z)  , 

(72) 

P  a2/2 (r2  +  F ( z ) )  . 

(73) 

Substitution  of  equations  (71)  through  (73)  into  equations 


(67)  through  (69)  yields: 
f'  2  (  2  )  -  2f  (Z)f  •  '  (z)  = 
2f  ( z ) f ' ( z )  = 


a2  +  vf '  '  ' 
a2F ' (z) /4  - 


(z)  , 

vf '  ' 


(z)  . 


(74) 

(75) 


54 


The  boundary  conditions  given  by  equations  (67a)  and  (68a) 
become : 


at  z  =  0,  f  =  0,  f  =  0,  and  F  =  0,  (74a) 

at  z  =  »,  f ’  =  a.  (75a) 

Using  the  similarity  transformation,  £=  (a/i/  )1/2Z  and 
f(z)  =  (a  V  ) 1/2  0 ( £  ) f  equation  (74)  reduces  to, 

0  1  '  '  +  200  ’  ’  -  0  ' 2  +  l  =  0,  (76) 

with  the  boundary  conditions  given  by, 

at  £  =  0 ,  0=0  and  0*  =  0,  (76a) 

at  £  -  oo,  0'  =1.  (76b) 

Equation  (76)  is  handily  solved  by  classical  Runge-Kutta 
methods.  Following  the  method  outlined  by  Hornbeck  {20}, 
let  xj  =  0  ,  X2  =  xj.',  and  X3  =  X2*.  By  doing 
so,  equations  (76)  ,  (76a) ,  and  (76b)  are  transformed  to  the 
following  matrix  system: 

X’  =  H(  |  ,X)  ,  (77) 

-  2x^x3  +  X2  -  1]T,  (78) 


H (  £  ,X)  =  [X2 
X(0)  =  XQ 


x3 

(0  0  K]T, 


(79) 


where  K  must  be  found  such  that  0*  ( « )  =  X2  (  »  )  =  1. 

The  algorithm  used  consists  of  the  following  steps: 


(1) 

Guess  a 

value  for  K, 

(2) 

Compute 

X*  =  Xn  +  A  t  /2H(  t 

(3) 

Compute 

M*  =  H(£n  +  A  t  /  2  , 

(4) 

Compute 

X**  =  Xn  +  A$/  2M*, 

(5) 

Compute 

M**  =  H  (  £  n  +  Af  /  2 

(6) 

Compute 

X***  =  xn  +  At  M**, 

(7) 

Compute 

M***  =  H(  $n  +  At 

(8) 

Compute 

M  =  1/6  tH<  t  w  *n) 

(9) 

Compute 

Xn+1  =  Xn  +  4*  M, 

(10)  Perform  steps  (2)  through  (9)  until  t  sufficiently 
approaches  infinity. 


(11)  If  X2  (  £  ao  )  is  sufficiently  close  to  1,  the 
calculation  is  complete, 


(12)  If  step  (11)  is  not  satisfied,  guess  a  new  value  for  K 
and  perform  the  above  steps  again. 


For  the  solution  of  equation  (76) ,  it  was  assumed  that  t  = 
10  is  sufficiently  close  to  infinity,  and  XQ  is  given  by 
equation  (79).  Using  At  =  0.01,  K  was  found  to  be  1.312. 
The  calculated  values  of  0,  <h' ,  and  '  are  plotted  in 


Fig.  6.  The  above  method  is  a  fourth  order  method,  which 
is  to  say  that  it  is  accurate  to  the  order  of  (A  $  )4. 

The  boundary  layer  thickness,  6,  is  traditionally 
defined  as  the  value  of  z  at  which  the  velocity  is  equal  to 
0.99  times  the  free  stream  value,  U.  In  terms  of  the 
quantities  plotted  in  Fig.  6,  =  u/U,  and  from  the 

Runge-Kutta  solution,  the  value  of  £  or  which  u/U  =  0.99  is 
about  1.98.  From  the  definition  of  £,  one  sees  that, 

(a/U  )1/2d  =  1.98.  (80) 

Knowing  the  stagnation  flow  constant,  a,  and  the 
kinematic  viscosity,  v,  the  boundary  layer  thickness  may 
be  calculated  from  equation  (80) .  The  kinematic  viscosity 
is  known  for  each  cell  from  the  dynamic  viscosity,  given  in 
Appendix  A,  and  the  fluid  density,  available  from  HYDRO. 

The  constant,  a,  however,  is  somewhat  troublesome  as  it  is 
not  actually  constant;  rather,  it  changes  with  time. 
Equation  (80)  is  a  steady-state  equation  which  indicates 
that  there  may  be  a  conflict  with  using  the  stagnation  flow 
solution.  It  is  felt  by  the  author  and  others  (3)  that  the 
problem  may  be  treated  as  pseudo-steady,  yielding 
qualitatively  good  results,  if  the  following  conditions  are 
met ; 

(1)  the  "constant,"  a,  settles  down  with  time, 

(2)  the  time  required  to  establish  the  boundary  layer 
is  small  compared  to  the  laser  pulse  time. 


CURVE  li  PHI 
CURVE  2:  PHI  1 
CURVE  3:  PHI ' 


Even  though  the  pseudo-steady  solution  is  an  approximation, 
it  should  identify  the  relative  role  of  convective  heat 
transfer  on  enhanced  thermal  coupling.  This  is  especially 
true  since  the  amount  of  numerical  data  available  for 
analysis  is  limited  to  about  one  third  of  the  pulse  time, 
and  many  results  are  projected  from  trends  during  latter 
times  of  the  data. 

The  time-dependent  constant,  a,  was  determined  from  the 
data  plotted  in  Fig.  5.  The  curves  in  this  plot  show  that, 
prior  to  the  LSD  front,  the  radial  velocity  component  is 
closely  approximated  by  equation  (66)  .  The  slope  of  the 
linear  portion  of  a  curve  for  time,  t,  is  a(t).  By  looking 
at  the  data  for  twelve  representative  times  (from  1.59 
microseconds  to  8.2  microseconds),  a(t)was  determined  and 
plotted  in  Fig.  7.  The  reason  only  twelve  times  were 
considered  is  that,  in  order  to  limit  the  amount  of  output, 
HYDRO  was  modified  to  output  data  at  every  20th  time  step 
(cycle),  and  due  to  limited  computer  resources,  only  twelve 
sets  of  output  were  possible.  Referring  to  Fig.  7,  curve  1 
is  the  computed  numerical  data  and,  though  the  data  goes 
only  to  8.2  microseconds ,  the  following  trend  was  noticed: 
a  ( 7 . 5)  s=  1/  2a  ( 7 . 0 )  ,  a(7.9)  =  l/2a(7.5),  and  a(8.2)» 
l/2a(7.9).  Based  on  this  trend,  it  was  projected  that 
a(t),  for  the  entire  pulse  time,  may  go  something  like 
curve  1.  Curve  2  is  a  forced  fit  to  curve  1,  chosen  so  as 
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to  approximate  the  general  trend  of  a(t).  The  equation  of 
curve  2  is : 

a ( t)  x  lO"6  =  0.76  exp(0. 35915  -  0.02339t  +  1.2118/t),  (81) 

for  times  less  than  18  microseconds,  and, 

a ( t)  x  10~6=  0.76,  (82) 

for  times  greater  than  18  micoseconds.  In  equation  (81),  t 
is  in  microseconds.  The  time-average  value  of  a(t)  was 
found  from  the  definition: 

a  =  1/rp  f  a (t) dt ,  (83) 

*o 

where  Tp  is  the  laser  pulse  time  and  a  is  the  time-averaged 
value  of  a(t).  From  equations  (81)  through  (83),  a  was 
found  to  be  about  0.95  x  10®  sec”l.  It  is  felt,  then, 
that  by  replacing  a,  in  the  steady-state  equations,  with  S, 
the  steady-state  equations  will  yield  qualitatively  good 
results  for  the  unsteady  flow.  Also,  if  the  projection  of 
a(t)  at  latter  times,  in  Fig.  7,  approximates  the  true 
behavior  (which  is  unknown  due  to  lack  of  data) ,  then 
condition  (1)  of  the  pseudo-steady  requirements  discussed 
earlier  is  satisfied,  in  that  a(t)  does  indeed  settle  down 


with  time. 


The  boundary  layer  thickness  is  given  by  equation 
(80) .  The  temperature  in  the  LSD  front,  at  the  target 
surface,  is  about  11,000  degrees  Kelvin.  From  Appendix  A, 
the  dynamic  viscosity  at  this  temperature  is  about  2.25  x 
10“3  gm/ (cm-sec) ,  and  the  density,  taken  from  HYDRO,  is 
about  0.95  x  10“4.  The  kinematic  viscosity,  V,  is  about 
23.7  (cm2/sec) ,  and,  from  equation  (80),  6  is  about  1.0  x 
10“2  cm.  The  time  required  to  establish  a  boundary  layer 
of  thickness,  <J,  may  be  approximated  by  the  relation  given 
by  Schlichting  (19)  as:  ts  6^/ 4V  .  it  is  interesting 
to  note  that,  according  to  equation  (80) ,  the  time  required 
to  establish  a  boundary  layer  in  3-D  stagnation  flow  is 
independent  of  flow  properties?  rather,  it  is  fully 
determined  by  the  stagnation  flow  constant  defined  by 
equation  (66) .  Using  this  relation,  the  time  required  to 
establish  the  boundary  layer  is  about  one  microsecond, 
which  is  about  four  percent  of  the  pulse  time,  25.73 
microseconds.  Thus,  condition  (2)  of  the  pseudo-steady 
requirements  is  satisfied,  in  that  the  time  required  to 
establish  the  boundary  layer  is  short  compared  to  the  laser 
pulse  time. 

With  the  above  fundamental  boundary  layer  concepts  in 
hand,  and  with  the  confidence  that  the  solution  may  be 
treated  as  pseudo-steady,  the  convection  heat  transfer  may 
now  be  calculated.  In  the  presence  of  a  LSD,  fluid 
temperatures  are  sufficiently  high  that  chemical 
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dissociation  occurs.  Consequently,  the  thermodynamic  and 
transport  properties  of  the  fluid  (air,  in  this  case)  vary 
greatly  across  the  boundary  layer.  For  these  types  of 
problems,  Eckert  £22}  has  recommended  use  of  the  reference 
enthalpy  technique  in  computing  heat  transfer.  Eckert 
reports  £22}  that  this  technique  yields  fairly  accurate 
results  (to  within  _+  20%  of  experimentally  measured 
values) ,  even  in  the  presence  of  compressibility,  fluid 
property  variation,  and  chemical  dissociation.  The 
reference  enthalphy  technique  is,  however,  a  steady-state 
technique.  Therefore,  before  using  it  to  calculate  the 
convection  heat  transfer  for  a  particular  flow  situation, 
one  must  verify  that  the  flow  may  be  approximated  as  being 
steady. 

The  radial  temperature  distribution  at  the  target 
surface,  as  given  by  HYDRO,  in  the  presence  of  the  LSD 
ignited  by  Hall's  9.3  J  pulse,  is  plotted  in  Fig.  8. 
Referring  to  this  plot,  one  sees  that  temperatures  on  the 
order  of  11,000  degrees  Kelvin  are  present?  therefore,  the 
reference  enthalpy  technique  was  chosen  to  perform  the 
convection  heat  transfer  calculations.  It  is  felt  that  use 
of  the  method  is  justified,  in  view  of  previous  arguments 
that  the  boundary  layer  flow  is  pseudo-steady  in  nature. 
Using  this  technique,  the  convection  heat  transfer 
coefficient  is  defined  by: 
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where  q  is  the  rate  of  heat  addition  to  the  target,  A  is 
the  area  over  which  the  heat  addition  is  taking  place,  h^ 
is  the  convection  heat  transfer  coefficient,  iw  is  the 
mass  specific  enthalpy  of  the  fluid  at  the  target  surface, 
and  ir  is  the  recovery  enthalpy  of  the  fluid  at  the  edge 
of  the  boundary  layer.  According  to  Eckert  {221,  the 
recovery  enthalpy,  in  the  neighborhood  of  a  3-D  stagnation 
point,  is  given  by. 


ir  =  ie  +  1/2  Pr3-/2  ue2. 


where  ie  is  the  mass  specific  enthalpy  at  the  edge  of  the 
boundary  layer,  ue  is  the  corresponding  fluid  velocity  and 
Pr,  the  Prandtl  number,  is  present  as  a  recovery  factor  for 
high  speed  flow.  Thus,  one  sees  that,  if  h^  is  known,  the 
amount  of  convection  heat  transfer  per  unit  area  is  easily 
determined  by  equation  (84),  with  ir  given  by  equation 
(85) .  The  problem,  therefore,  consists  primarily  of  finding 
h^.  According  to  Eckert  {22},  if  the  external  flow  field 
is  that  of  3-D  stagnation  flow,  h^  may  be  found  from  the 
following  expression: 


Nui  =  0.763 (Pri) °*4 (Rei) °*5, 


where  Nu^  is  the  Nusselt  number,  Pr^  is  the  Prandtl 
number,  and  Re^  is  the  Reynolds  number.  In  equations 
(84)  and  (86)  ,  the  subscript,  i,  denotes  that  the 
properties  must  be  evaluated  at  the  reference  enthalpy 
temperature  (Tref)  which  is  determined  from  the  following 
relation: 


=  ie  +0 . 5  (iw  ”  ie)  +  0 . 22  ( ir  i0 )  ,  (87) 

where  iref  is  the  reference  enthalpy.  The  nondimensional 
numbers  appearing  in  equation  (86)  are  defined  as  follows: 

Nui  =  Cph^r/k^,  (88) 

i 

Prj_  =  //iCp/kj_,  (89) 

i 

Rei  =  uer/ v  ir  (90) 

where  r  is  the  radius  from  the  axis  of  symmetry.  Note  that 
the  definition  of  the  Nusselt  number  in  equation  (88) 
differs  from  the  traditional  definition  by  the  factor, 

Cp..  This  is  due  to  the  fact  that  hj_  is  defined  in 
terms  of  enthalpy  differences  rather  than  temperature 
differences.  From  equation  (66),  ue  =  ar,  where  a  has 
been  replaced  by  the  time-average  value  of  a(t),  found 
previously.  With  this  result,  equations  (85)  and  (90) 


become : 


1/2  Pr 1/2 (ar) 2, 


(91) 


ir  =  ie+ 

Rei  =  ar2/ v  (92) 

Susbstitution  of  equations  (88) ,  (89) ,  and  (92)  into 
equation  (86)  yields: 

hi  =  0.763(5)1/2  Pil/2  /ii-0.1  (ki/c  v0.6f  (93) 

*1 

where  a  =  0.95  xlO6 7 8  sec“l. 

The  algorithm  for  computing  the  convection  heat 
transfer  is  now  complete  and  consists  of  the  following 
steps : 

(1)  Obtain  T  and  pressure  from  HYDRO, 

(2)  Compute  ie  from  the  curve  fit  given  in  Appendix  A, 

(3)  Compute  ir  from  equation  (91), 

(4)  Compute  iref  from  equation  (87) , 

(5)  Compute  Tref  by  iteration,  such  that  the  curve 
fit  yields  iref, 

(6)  Compute  p  i  and  Cp.  at  Tref,  using  the 

curve  fits  given  in  Appendix  A, 

(7)  Compute  hi  from  equation  (93) , 

(8)  Compute  q/A  from  equation  (84)  . 


The  above  algorithm  is  used  to  compute  the  heat  transfer 
out  to  the  radius  at  which  the  temperature  is  approximately 
equal  to  the  ambient  temperature,  taken  to  be  the  standard 
temperature  of  0°C.  The  target  temperature  is  taken  to  be 
the  same  standard  temperature.  The  last  assumption  is  not 
correct  since,  if  target  vaporization  is  taking  place,  the 
surface  temperature  in  these  regions  is  essentially  at  the 
vaporization  temperature  appropriate  for  the  local 
pressure;  this  will  always  be  much  higher  than  O^C. 
Therefore,  the  heat  transfer  obtained  with  the  assumption 
should  be  higher  than  the  true  value.  Also,  it  is 
important  to  note  that  in  front  of  the  LSD,  where  flow 
velocities  are  small,  the  actual  values  of  the  local  edge 
velocity  given  by  HYDRO  are  used  in  computing  the  heat 
transfer.  This  correction  to  the  above  algorithm  is 
necessary  since  ue,  as  defined  by  equation  (66),  is  valid 
only  to  the  LSD  front. 

Due  to  limited  computer  resources,  the  heat  transfer 
calculations  are  done  by  a  separate  program,  for  which  the 
input  data  is  the  output  from  HYDRO.  This  program  is 
listed  in  Appendix  C. 

In  summary,  it  is  necessary  to  assess  the  flow 
situation  and  to  predict  a  few  fundamental  boundary  layer 
quantities  before  one  chooses  a  method  by  which  to  compute 
the  amount  of  convection  heat  transfer  taking  place.  The 
determination  of  whether  the  boundary  layer  flow  may  or  may 
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not  be  approximated  as  steady-state  is  of  particular 
importance,  as  most  known  solutions  are  for  steady  flow. 
Also,  it  is  important  to  predict  the  effect  of 
compressibility  on  the  boundary  layer,  as  flow  properties 
may  change  drastically  across  the  boundary  layer.  In  this 
section,  it  was  deduced  from  the  hydrodynamic  data  that  the 
external  flow  field  for  the  problem  under  study  closely 
resembles  that  of  three-dimensional  axisymmetric  stagnation 
flow.  The  boundary  layer  thickness  and  the  time  required 
to  establish  the  boundary  layer  were  estimated  from  the 
compressible,  steady  flow  boundary  layer  solutions.  From 
these  considerations,  it  was  concluded  that  the  boundary 
layer  flow  may  be  treated  as  pseudo-steady.  High  fluid 
temperatures  and  the  possible  effects  of  compressibility 
and  property  variations  on  the  heat  transfer  solution  led 
to  the  choice  of  Eckert's  reference  enthalpy  technique  as 
the  method  used  to  compute  the  heat  transfer,  since  this 
technique  is  known  to  give  fairly  good  results  in  the 
presence  of  such  factors.  Finally,  the  algorithm  used  to 
compute  the  heat  transfer  was  presented.  The  results  of 
these  calculations  are  presented  in  Section  VII. 
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V.  Vapor  Condensation 


As  stated  in  the  introduction,  the  second  object  of 
this  study  is  to  examine  the  role  of  vapor  condensation  on 
enhanced  coupling.  In  this  section,  the  method  used  to 
calculate  the  amount  of  thermal  energy  imparted  to  the 
target  by  condensing  target  vapor  is  discussed  in  detail. 
The  motivation  for  performing  these  calculations  was  stated 
in  the  introduction  and  is  not  repeated  here.  It  is  noted, 
however,  that  the  presence  of  a  high-pressure  LSD  is 
required  for  vapor  condensation  to  take  place  because  it  is 
the  high  pressures  associated  with  the  LSD  that  force  the 
evolved  target  vapor  back  onto  the  surface.  In  this 
section,  therefore,  the  presence  of  a  LSD  is  inherently 
assumed. 

When  a  mixture  of  air  and  target  vapor  passes  over  a 
surface,  the  vapor  species  is  deposited  so  as  to  remove 
energy,  mass,  and  momentum  from  the  flow.  According  to 
Jackson  [24,  25},  these  quantities  are  determined  from 
simple  Maxwell-Boltzmann  statistics.  The  Maxwellian 
velocity  distribution  function  is  given  by  Vincenti  and 
Kruger  {15}  to  be: 


where  m  is  the  mass  of  the  vapor  particles,  T  is  the 
temperature  of  the  vapor  particles,  is  the  radial 
velocity  component,  C2  is  the  velocity,  component  in  the 
angular  direction,  and  C3  is  the  axial  velocity 
component.  Equation  (94)  gives  the  probability  that  a 
given  molecule  chosen  at  random  will  have  velocity  in  the 
range  Ci ,  C2>  C3  to  +  dC]_,  C2  +  dC2,  C3 
+  dC3*  It  is  assumed  that,  for  this  problem,  vapor 
particles  immediately  above  the  target  surface  impinge  upon 
the  surface  with  the  instantaneous  velocity,  C3,  relative 
to  the  surface.  Furthermore,  it  is  assumed  that  all  vapor 
particles  which  come  into  contact  with  the  surface  will  be 
attached  to  it  where  they  will  condense  and  return  to  the 
target  the  energy  required  to  vaporize  them  as  well  as  the 
thermal  kinetic  energy  gained  while  in  the  vapor  state. 

The  last  assumption  needs  verification  and,  according  to 
Levine  and  Gyftopoulos  (26),  virtually  every  particle  which 
comes  into  contact  with  the  surface  does  attach  to  the 
surface,  where  it  condenses. 

With  the  above  assumptions,  the  flux  of  quantity,  Fj, 
into  the  target  surface  may  be  calculated  from  the 
following  relations: 


where  n  is  the  particle  number  density,  Fj  is  the 
quantity  imparted  to  the  surface  by  vapor  collisions, 
f (CjJ  is  given  by  equation  (94),  and  the  axial  velocity 
goes  from  zero  to  - «  since  it  is  zero  at  the  surface  and 
the  impingement  is  in  the  negative  axial  direction.  The 
Fj  and  Fj  in  equation  (95)  are  given  by: 


=  mass  flux,  m, 

=  energy  flux,  q, 

=  normal  momentum  flux,  M3, 

=  radial  momentum  flux,  M4, 
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-  kinetic  energy,  mC2/2, 
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The  appearance  of  the  latent  heat  of  vaporization,  Lv,  in 
equation  (97)  is  due  to  the  fact  that  condensing  particles 
must  return  the  energy  originally  required  to  vaporize 
them.  The  mass  flux  of  particles  is  calculated  by 
inserting  equation  (94)  into  equation  (96), 
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=  n  y*(m/2  WkT)  3/2  mC3exp (-mC3/ 2kT) dC3  l\, 
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where  1^  =  exp (-mx2/2kT) dx* .  1^  is  given  by 

Vincenti  and  Kruger  [153  to  be  ( 2  It  kT/m)  1/ 2 .  Equation 


(96a)  becomes, 
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m  =  mn 


(2  7rkT/m)  (m/2  n  kT)  3/2  J  c3exp  (-mC^/ 2kT)  dC3  . 


(96b) 


Performing  the  integration,  equation  (96b)  becomes. 


m  =  mn  (  2  ft  kT/m)  (m/  2  7T  kT)  3/ 2kT/m. 


(96c 


The  term,  mn,  is  defined  to  be  the  vapor  density; thus ,  the 
mass  flux  is  given  by: 

m  =  Pv(kT/2  7Tm)  3-/2 .  (100 
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In  this  section,  x  and  y  are  dummy  integration  variables 


The  energy  flux  imparted  to  the  surface  by  condensing  vapor 
particles  is  determined  by  inserting  equation  (94)  into 
equation  (97) , 
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ALV  III  C3exp [B (xz  +y^  +C3 ) ] dydxdC3 ,  (97a) 
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where  A  =  mn  (m/  2  n  kT)  3/  2  =  (m/ 2  ?r  kT)  3/  2  and  B  =  -m/2kT. 

Denoting  the  first  triple  integral  in  equation  (97a)  by 
l2' 
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Kruger  {15}  to  be  (2  7T  kT/m)  1/ 2 .  Equation  (97b)  becomes, 
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The  integral  in  equation  (97c)  is  given  {15}  to  be 
1/2B2  .  Thus,  with  the  insertion  of  A  and  B,  equation 

(97c)  becomes, 


1 2  =  A/  (4  *r2)  (2  *r  kT/m)  3. 


( 9  7d) 


Denoting  the  second  triple  integral  in  equation  (97a)  by 
14- 
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where  I5  =  x2exp  Bx2dx,  given  (153  to  be  1/(2  ir)(2  7TkT/m) 
Evaluating  the  integral  in  equation  (97e)  and  inserting  the 
above  values  for  I3,  I5,  A,  and  B,  equation  (97) 
becomes , 


I4  =  A/<4  n2)  (2rkT/m)3.  (97f) 

The  third  triple  integral  in  equation  (97a)  is  simply 
mLv;  therefore,  the  energy  flux  is  determined  to  be, 

q  =  m ( 2kT/m  +  Ly) .  (101) 

An  approach  similar  to  that  taken  above  shows  that  M3  and 
M4  are,  respectively,  pv! 2  (kT/m)  and  zero.  The 
derivation  is  not  included  here,  however,  since  the  amount 
of  thermal  coupling  is  determined  completely  by  equations 
(1 00)  and  (101). 

In  summary,  it  is  felt  that  vapor  particles  passing 
over  the  target  surface  will  become  attached  to  it;  thereby 
imparting  the  flux  of  mass,  energy,  and  momentum  to  the 


target.  The  energy  contained  by  the  particles  is  believed 
to  be  greater  at  the  time  of  condensation  than  at  the  time 
of  evolution  due  to  the  fact  that  they  gain  kinetic  energy 
while  in  the  vapor  state.  In  this  section,  simple 
Maxwel 1-Bol tzmann  kinetic  theory  is  used  to  evaluate  the 
amount  of  mass,  energy,  and  momentum  imparted  to  the 
surface  by  the  condensing  particles.  The  fluid  properties 
needed  for  the  calculations  are  available  from  HYDRO  and 
the  results  of  these  calculations  are  presented  in  Section 


VI .  Radiation  Heat  Transfer 
The  third  and  final  object  of  this  study  is  to  estimate 
the  effect  of  plasma  radiation  on  enhanced  coupling.  When 
a  LSD  is  ignited  at  the  target  surface,  it  propagates  up 
the  laser  beam  and,  due  to  its  high  absorptivity,  it 
absorbs  practically  all  of  the  incident  laser  radiation. 

It  is  assumed  that  some  or  all  of  this  absorbed  radiation 
will  be  radiated  away  from  the  plasma  by  various 
mechanisms.  This  process  is  referred  to  as  plasma 
reradiation.  The  problem  is  to  describe  which  processes 
are  responsible  for  plasma  reradiation  and  to  estimate  the 
amount  of  energy  imparted  to  the  target  by  these  processes. 

According  to  Dawson  (27) ,  there  are  two  regimes  of 
plasma  reradiation:  bremsstrahlung  radiation  and  blackbody 
radiation.  At  a  given  point  in  space  and  time,  the  type  of 
radiation  present  is  determined  by  the  absorption  length  of 
the  plasma.  The  absorption  length,  1,  is  by  definition  the 
inverse  of  the  free-free  absorption  coefficient  given  by 
equation  (37)  of  Section  III.  According  to  Dawson,  if  1  is 
less  than  the  plasma  radius  (or  thickness  for  a  planar 
plasma)  the  radiation  is  blackbody  and,  if  1  is  greater 
than  the  plasma  radius  (or  thickness)  the  radiation  is 
bremsstrahlung. 

Cjl 


For  bremsstrahlung  radiation,  the  volumetric  rate  of 
radiation  is  given  by  Dawson  (27)  to  be, 

€=  4.86  x  10“31Zn|  T3-/2  (102) 

where  f  is  the  volumetric  rate  in  watts  per  cubic 
centimeter,  Z  is  the  degree  of  ionization,  ne  is  the 
electron  number  density,  and  T  is  the  temperature  in 
kiloelectron  volts.  The  maximum  temperature  in  the 
numerical  mesh,  as  given  by  HYDRO,  is  only  about  0.002 
kiloelectron  volts  for  the  problem  under  study;  therefore, 
the  degree  of  ionization  is  assumed  to  be  one.  Actually, 
this  assumption  is  implicit  in  that  HYDRO  uses  the  Saha 
equation  for  singly  ionized  species  in  the  derivation  of 
the  electron  number  density.  With  this  assumption, 
equation  (102)  becomes, 

€  =  1.427  x  10“34n|T1/2.  (103) 

where  T  is  now  in  degrees  Kelvin.  Inserting  equation  (37) 
into  equation  (103)  yields: 


€=  5.715  x  10“6KffT2. 
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Assuming  that  the  LSD  front  is  planar,  its  volume  is  simply 
its  thickness  times  its  cross-sectional  area.  From  HYDRO, 
it  appears  that  the  thickness  of  the  LSD  is  on  the  order  of 
0.1  centimeters.  Actually,  the  thickness  is  expected  to  be 
less  than  that  deduced  from  HYDRO  due  to  the  fact  that  the 
artificial  diffusion  has  the  effect  of  spreading  the  front 
such  that  it  appears  to  be  thicker  than  it  really  is.  It 
is  felt,  therefore,  that  using  the  thickness  predicted  by 
HYDRO  yields  the  maximum  radiation  possible.  With  the 
above  assumption,  the  heat  flux  leaving  the  plasma  due  to 
bremsstrahlung  radiation  is, 

q  =  5.715  X  10_7kffT2,  (105) 

where  q  is  the  heat  flux  in  watts  per  square  centimeter.  It 
is  not  expected  that  all  of  the  radiation  which  leaves  the 
plasma  will  reach  the  target  surface;  thus,  the  inclusion  of 
a  radiation  shape  factor,  as  discussed  by  Holman  (21), 
should  be  made  in  equation  (105)  .  Furthermore,  it  is  not 
expected  that  all  of  the  radiation  incident  on  the  target 
surface  will  be  absorbed;  therefore,  the  inclusion  of  a 
surface  absorptivity  should  be  made  in  equation  (105) .  As 
discussed  in  Section  III,  the  absorptivity  of  the  surface 
may  be  assumed  to  be  unity  in  very  early  times  due  to  the 
mechanism  of  unipolar  arcing.  However,  once  a  LSD  is 
formed,  the  absorptivity  of  the  surface  rapidly  decreases. 


Based  on  these  ideas,  it  is  felt  by  the  author  that  the 
heat  flux  absorbed  by  the  surface  due  to  bremsstrahlung 
radiation  is  greatly  overestimated  by  equation  (105)  .  In 
fact,  it  is  suspected  that  the  predicted  value  may  be  an 
order  of  magnitude  greater  than  the  true  value.  Thus, 
equation  (105)  should  yield  the  maximum  possible  value  of 
thermal  coupling  due  to  bremsstrahlung  radiation  from  the 
LSD  front.  It  is  noted  that  the  above  model  neglects 
reradiation  from  bound  electrons  and  from  recombination 
but,  according  to  Dawson  (27),  these  processes  are 
unimportant  for  plasma  temperatures  below  the  order  of  one 
hundred  electron  volts.  Since  the  maximum  temperature  for 
this  problem  is  on  the  order  of  two  electron  volts,  it  is 
felt  that  these  processes  do  not  contribute  to  the  plasma 
radiation. 

According  to  Dawson,  (27},  blackbody  radiation  is 
negligible  for  plasma  temperatures  below  the  order  of 
several  hundred  electron  volts;  therefore,  it  is  not 
addressed  in  this  study. 

In  summary,  the  type  of  plasma  reradiation  present  is 
determined  by  the  absorption  length  of  the  LSD  front.  The 
model  presented  by  Dawson  (27}  was  used  to  estimate  the 
heat  flux  absorbed  by  the  target  due  to  bremsstrahlung 
radiation  from  the  plasma,  given  by  equation  (105).  It  is 
felt  that  this  equation  will  greatly  overestimate  the 


amount  of  radiation  absorbed,  due  to  the  absence  of  a 
surface  absorptivity  and  a  radiation  shape  factor  in  the 
equation.  Blackbody  radiation  is,  according  to  Dawson 
[27},  negligible  for  the  problem  under  study.  The  results 
of  the  bremsstrahlung  radiation  estimates  are  presented  in 


Section  VII. 


VII.  Results 


In  this  section,  the  results  of  the  calculations 
discussed  in  the  previous  sections  are  presented.  It  is  of 
utmost  importance  to  note  that  the  results  of  this  study 
are  based  on  a  limited  amount  of  numerical  data,  as  stated 
in  the  introdution.  Nevertheless,  the  author  feels  that 
the  results  are  qualitatively  correct  and  that  the  relative 
roles  of  convection  heat  transfer,  vapor  condensation,  and 
radiation  heat  transfer  on  enhanced  coupling  are  accurately 
predicted  from  these  results. 

The  amount  of  convection  heat  transfer  taking  place  at 
the  target  surface  under  the  presence  of  a  LSD  was 
calculated  from  the  algorithm  given  in  Section  IV.  The 
resulting  heat  flux  into  the  surface  is  plotted,  as  a 
function  of  radius  from  the  axis  of  symmetry  with  time  as  a 
parameter,  in  Fig.  9.  Referring  to  Fig.  9,  one  sees  that 
at  early  times,  prior  to  the  ignition  of  a  LSD,  the 
convection  heat  transfer  is  quite  small.  This  result  is 
consistent  with  the  theory  that  it  is  the  presence  of  the 
LSD  which  facilitates  enhanced  coupling.  After  the  LSD  is 
ignited  (at  about  5  microseconds,  as  discussed  later),  one 
sees  from  Fig.  9  that  the  convection  heat  transfer 
increases  by  a  factor  of  about  7.  As  time  progresses, 
however,  one  sees  that  the  convection  heat  transfer  appears 
to  converge  to  a  radial  distribution  with  no  dependence  on 
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time.  The  time  for  this  convergence  is  at  about  7 
microseconds.  Based  on  these  observations,  it  is  projected 
that  the  convection  heat  transfer  may  be  approximated  by  a 
single  radial  distribution  for  all  times  up  to  8.2 
microseconds,  which  is  taken  to  be  about  10  microseconds . 
This  distribution  is  plotted  in  Fig.  10.  From  Fig.  10,  one 
sees  that  the  projected  heat  flux  overestimates  the 
early-time  curves  of  Fig.  9,  underestimates  the  middle-time 
curves  of  Fig.  9,  and  approximates  the  late-time  curves  of 
Fig.  9.  To  obtain  the  energy  imparted  to  the  surface  per 
unit  area,  the  heat  flux  curve  of  Fig.  10  was  integrated 
over  time,  up  to  10  microseconds.  Since  the  curve  is 
independent  of  time,  the  integration  amounts  to  simply 
multiplying  the  curve  by  10  microseconds.  Looking  ahead, 
the  result  is  plotted  as  curve  2  in  Fig.  13  (p.  90) ,  from 
which  one  sees  that:  (a)  the  amount  of  energy  absorbed  by 
the  target  due  to  convection  heat  transfer  in  10 
microseconds  is  much  less  than  that  measured  by  Hall  (curve 
1) ,  and  (b)  the  radial  distribution  of  the  energy  absorbed 
due  to  convection  heat  transfer  is  not  the  same  as  for  curve 
1  (Fig.  13  p.  90) .  The  thermal  coupling  coefficient,  at  a 
given  time,  is  simply  the  thermal  energy  residing  in  the 
target  at  that  time  divided  by  the  energy  in  the  laser  beam 
at  that  time.  Therefore,  the  thermal  coupling  due  to 
convection  heat  transfer,  at  10  microseconds,  is  curve  2 
(Fig.  13  p.  90)  integrated  over  the  area,  divided  by  the 


beam  energy  at  10  microseconds.  Curve  2  was  (Fig.  13  p. 
90)  approximated  by  the  following  curve  fit, 

E/A*  0 . 444exp [a (r) ] ,  (106 

where  E/A  is  the  energy  per  unit  area  and, 

a(r)=30.7r2  -  641r3  +  2231r4  -  2405r5.  (107) 

The  total  energy  absorbed  by  the  target  is, 

0.4 

Et=s  2  IT  J  E/A  rdr .  (108 

0 

Equation  (108)  was  integrated  numerically  by  the  trapezoid 
rule  using  Ax  =  0.0004,  with  the  result  that 
Et  =  0.0558  J.  From  Fig.  3,  the  beam  energy  at  10 
microseconds  is  about  5  J;  therefore,  the  thermal  coupling 
due  to  convection  heat  transfer  at  10  microseconds  is  only 
about  0.011. 

The  energy  flux  carried  into  the  target  by  condensing 
vapor  particles  was  calculated  from  equation  (101)  of 
Section  V.  This  energy  flux  is  plotted,  as  a  function  of 
radius  from  the  axis  of  symmetry  with  time  as  a  parameter, 
in  Fig.  11.  Note  that  the  curves  in  Fig.  11  are  plotted 
only  for  times  after  the  ignition  of  the  LSD  (about  5 
micoseconds)  since  equation  (101)  assumes  that  all  vapor 
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particles  which  come  into  contact  with  the  surface  condense 
and  this  is  not  true  prior  to  formation  of  the  LSD,  due  to 
the  fact  that  surface  vaporization  takes  place  until  the 
LSD  forms.  From  Fig.  11,  one  sees  that  the  energy  flux 
into  the  surface  due  to  vapor  condensation  decreases  with 
time  and,  as  in  the  case  of  convection  heat  transfer,  the 
flux  appears  to  converge  to  a  radial  distribution  with  no 
time  dependence  by  about  7  microseconds.  It  is  projected, 
then,  that  the  energy  flux  may  be  approximated  by  a  single 
radial  distribution,  independent  of  time,  for  all  times  up 
to  about  10  microseconds.  This  distribution  is  plotted  in 
Fig.  12,  from  which  one  sees  that  the  projected  energy  flux 
underestimates  the  early-time  curves  of  Fig.  11, 
approximates  the  middle-time  curves  of  Fig.  11,  and 
overestimates  the  late-time  curves  of  Fig.  11.  To  obtain 
the  energy  imparted  to  the  target  per  unit  area,  the  energy 
flux  of  Fig.  12  was  integrat -.d  over  time  up  to  10 
microseconds.  As  before,  the  integration  amounts  to  simply 
multiplying  the  curve  in  Fig.  12  by  10  microseconds.  The 
result  is  plotted  as  curve  3  in  Fig.  13,  from  which  one 
sees  that:  (a)  the  energy  imparted  to  the  target  by  vapor 
condensation  is  much  greater  than  the  energy  imparted  to 
the  target  by  convection  heat  transfer,  (curve  2)  and  (b) 
the  measured  radial  distribution  of  the  energy  imparted  to 
the  target  per  unit  area  (curve  1)  is  more  closely  matched 
by  the  condensation  curve  than  by  the  heat  transfer  curve. 
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In  order  to  find  the  thermal  coupling  due  to  vapor 
condensation,  curve  3  was  fit  to  the  following  expression: 

E/A  =  1.2exp(-4.76r2) .  (109) 

Integrating  equation  (112)  over  the  area,  the  energy  residing 
in  the  target  due  to  vapor  condensation  is, 

0.7 

Et  s:  27 T  J E/A  rdr.  (110) 

0 

Equation  (110)  gives  Et  to  be  about  0.715J.  The  energy 
in  the  beam  at  10  microseconds  is  about  5  J;  therefore,  the 
thermal  coupling  due  to  vapor  condensation  at  10 
microseconds  is  about  0.143.  Curve  4  of  Fig.  13  shows  the 
additive  effect  of  convection  heat  transfer  and  vapor 
condensation  on  the  energy  absorbed  by  the  target. 

The  radiation  heat  transfer  to  the  target  from  the  LSD 
front  was  estimated  from  equation  (105)  of  Section  VI.  If 
the  absorption  coefficient  defined  by  equation  (37)  of 
Section  III  is  greater  than  10,  the  radiation  is  blackbody 
and  is  neglected  according  to  Dawson's  127)  recommendation 
discussed  in  Section  III.  Otherwise,  the  radiation  is 
bremssrahlung  and  equation  (105)  applies.  From  the 
numerical  data  given  by  HYDRO,  it  was  observed  that,  prior 
to  ignition  of  the  LSD,  the  radiation  is  bremsstrahlung, 


and  afterwards,  the  radiation  is  blackbody.  The 
calculation  of  the  heat  flux  into  the  target  due  to 
bremsstrahlung  radiation  is  summarized  in  Table  I. 


Table  I 

Estimate  of  Bremsstrahlung  Plasma  Reradiation 


t  x  10“6  sec 

kff  cm-1 

T  (deg  K) 

q  W-cm-2 

1.59 

2.31 

2942 

11.43 

2.71 

2.24 

3747 

17.97 

3.73 

1.95 

3767 

15.81 

From  Table  I,  one  sees  that  the  maximum  heat  flux  into 
the  target  due  to  bremsstrahlung  radiation  is  about  18 
watts  per  square  centimeter.  For  reasons  discussed  in 
Section  VI,  it  is  felt  that  this  number  greatly  over¬ 
estimates  the  true  radiation  heat  flux  into  the  target. 
Assuming  that  the  radiation  heat  flux  into  the  target  is 
the  maximum  value  given  above  and  that  it  is  constant 
through  the  entire  laser  pulse,  one  obtains  that  the 
energy  absorbed  by  the  target  per  unit  area  is  only  about 
4.68  x  10“4  J/cm2.  Furthermore,  if  one  assumes  that 
this  value  is  constant  over  the  entire  laser  spot  radius, 
one  obtains  that  the  total  energy  transferred  to  the 
target  by  radiation  is  about  5.9  x  10-5  J.  The  beam 
energy  at  the  end  of  the  pulse  is  9.3  J;  therefore,  the 

92 


JL*.»  'A— A  ‘A  ■  /  '-.S  1  V 


PARAMETER*  TIME  IN 
MICROSECONOS 


14:  DEPTH  OF  LASER- I NDUCED  VAPORIZATION 


thermal  coupling  due  to  bremsstrahlung  radiation  from  the 
LSD  front  is  only  6  x  10“5  percent.  This  number  is 
extremely  small  compared  to  the  coupling  via  convection 
heat  transfer  and  vapor  condensation.  Note  also  that  this 
extremely  small  number  is  felt  to  greatly  overestimate  the 
true  value.  It  appears  that  plasma  radiation  is  an 
insignificant  coupling  mechanism  for  the  problem  under 
study.  This  is  consistent  with  results  given  by  Jumper 
[28} ,  in  which  he  accurately  predicted  the  propagation 
speed  of  a  LSD  using  a  global  energy  balance  and  completely 
ignoring  radiation  away  from  the  LSD. 

It  was  stated  above  that  the  time  of  LSD  ignition  is 
about  5  microseconds.  Referring  to  Fig.  14,  where  the 
depth  of  laser- induced  target  vaporization  is  plotted,  one 
sees  that  by  about  5  microseconds,  laser  radiation  has 
practically  ceased  reaching  the  target.  The  reason  for  the 
rapid  decrease  in  vaporization  is  that  the 
surface-shielding  LSD  has  formed,  detached  from  the 
surface,  and  is  absorbing  practically  all  of  the  incident 
radiation.  Thus,  it  is  seen  that  the  LSD  is  formed  at 
about  5  microseconds. 

In  summary,  the  energy  imparted  to  the  target  per  unit 
area  by  convection  heat  transfer  was  calculated  with  the 
methods  of  Section  IV.  The  results  of  these  calculations 
are  presented  in  Fig.  9  and  Fig.  10.  The  energy  imparted 
to  the  target  per  unit  area  by  vapor  condensation  was 


Of. 


calculated  with  the  methods  of  Section  V.  The  results  of 
these  calculations  are  presented  in  Fig.  11  and  Fig.  12. 

The  relative  magnitudes  (and  radial  distributions)  of  these 
calculations  are  compared  both  with  each  other  and  with 
Halls  experimental  profile  in  Fig.  13.  Finally,  the 
radiation  heat  flux  into  the  target  was  estimated  using 
Dawson's  recommendations  discussed  in  Section  VI.  The 
results  of  these  estimations  are  shown  in  Table  I. 
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VIII.  Conclusions 


The  conclusions  of  this  study  were  made  based  on  the 
results  discussed  in  Section  VII.  From  the  results  shown 
in  Fig.  9  and  Fig.  10,  it  was  concluded  that  convection 
heat  transfer  is  only  a  minor  contributor  to  enhanced 
coupling  observed  by  Hall  for  a  9.3  J  laser  pulse.  First 
of  a.’l,  convection  heat  transfer  does  not  give  the  measured 
magnitude  of  energy  absorbed  by  the  target;  furthermore,  it 
does  not  give  the  correct  radial  distribution  of  the 
absorbed  energy.  From  the  results  shown  in  Fig.  11  and 
Fig.  12,  it  was  concluded  that  vapor  condensation  could 
essentially  account  for  most  of  the  enhanced  coupling, 
since  the  magnitude  of  energy  "returned"  to  the  target  is 
not  only  higher  than  that  for  the  heat  transfer,  but 
furthermore,  the  radial  distribution  of  the  absorbed  energy 
more  closely  resembles  Hall's  measured  profile,  as  shown  m 
Fig.  13.  The  above  results  are  based  on  numerical  data 
spanning  10  microseconds  of  a  25  microsecond  pulse. 

However  the  author  feels  that  the  projected  curves  shown  in 
Fig.  10  and  Fig.  12  will  not  increase  much,  if  any,  for 
times  greater  than  10  microseconds  due  to  the  fact  that  the 
incident  laser  radiation  falls  to  zero  fairly  rapidly  past 
this  time  (see  Fig.  2).  Therefore,  the  above  conclusions 
are  felt  to  be  justified  for  the  entire  pulse.  Finally,  it 
was  concluded  that  plasma  reradiation  plays  practically  no 
role  in  enhanced  thermal  coupling. 


The  above  conclusions  are  for  a  particular  problem, 
that  of  Hall's  9-3  J  pulse.  It  is  recommended  that  further 
studies  in  this  area  be  directed  towards  examining  the 
effects  of  laser  beam  and  target  material  parameters  on 
enhanced  thermal  coupling. 


Appendix  A  ^ 

Thermodynamic  and  Transport  Properties  of 
High  Temperature  Air 

The  curve  fits  presented  in  this  appendix  were  derived 
from  tabulated  data  given  by  reference  [29},  using  the 
method  of  least  squares  given  in  reference  [13}.  All  of 
the  fits  are  polynomials  and,  though  they  are  somewhat 
tedious,  they  accurately  match  Hansen's  [29}  tabulated 
data.  Curve  fits  were  obtained  for  the  following 
properties:  compressibility  factor  (Z) ,  mass  specific 
enthalpy  (i)  in  calories  per  gram,  mass  specific  heat  at 
constant  pressure  (cp)  in  calories  per  gram  per  degree 
Kelvin,  dynamic  viscosity  ( /i  )  in  grams  per  centimeter  per 
second,  thermal  conductivity  (k)  in  calories  per  centimeter 
per  second  per  degree  Kelvin,  and  Prandtl  number.  The 
curves  are  given  as  functions  of  temperature  (times  1/1000 
degrees  Kelvin)  with  pressure  (in  atmospheres)  as  a 
parameter.  With  the  pressure  and  temperature  given  by 
HYDRO,  the  properties  are  found  in  the  following  manner: 
if  the  pressure  is  less  than  one  atmosphere,  the  one 
atmosphere  value  is  used;  if  the  pressure  is  greater  than 
ten  atmospheres,  the  ten  atmosphere  value  is  used  and;  if 
the  pressure  is  between  one  and  ten  atmospheres, 
logarithmic  (to  base  10)  interpolation  is  used  to  find  the 


properties . 


The  actual  curve  fits  are  not  included  here  due  to  the 
fact  that  they  are  quite  elementary  in  nature  and  very 
bulky,  as  some  of  them  are  fourth  and  fifth  order 
polynomials . 


IE-3  CRL/GM 


Cm  FIG.  fl.2:  MASS  SPECIFIC  ENTHALPY  FOR  AIR 
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FIG.  A. 3:  SPECIFIC  HEAT  AT  CONSTANT  PRESSURE  ( 


1000K  CflL/( CM-SEC-K ) 
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Variable 

Definition 

Units 

A(ANE) 

Constant  used  in  finding  electron 
number  density 

- 

A(ENV) 

Constant  used  in  defining  spatial 
variation  of  a  Gaussian  beam 

1/cm2 

ABC 

Total  absorption  coefficient 

1/cm 

ABCC 

"Cold"  absorption  coefficient 

1/cm 

ABCF 

Frge-free  absorption  coefficient 

1/cm 

ABCM 

Maximum  absorption  coefficient  in 
mesh 

1/cm 

ABCO 

Constant  used  in  finding  ABCC 

cm2/g 

ABF 

Constant  used  in  finding  ABCF 

- 

AJ 

Area  of  Jth  radial  cell 

cm2 

ALF 

Factor  used  in  forward/backward 
differencing 

- 

AN 

Array  of  hydrodynamic  quantities 

- 

ANE 

Electron  number  density 

1/cm3 

ANEL 

Electron  number  density 

1/cm3 

B(MHD) 

Factor  used  in  forward/backward 
differencing 

- 

B (ANE) 

Constant  used  in  finding  electron 
number  density 

- 

B(ENV) 

Exponent  in  exp(B)  for  a  Gaussian 
beam 

- 

C 

Cosine  0  ,  used  in  finding  electron 
number  density 

- 

CA 

Cos-1  0  ,  as  defined  above 

- 

CB 

Constant  used  in  finding  electron 
number  density 

- 

CF1 


Space-interval  pressure  correction 
factor 


Variable 


Definition 


Units 


CF2 

Time-interval  pressure  correction 
factor 

CHEW 

Depth  of  laser-induced  target 
vaporization 

cm 

CTV 

Product  of  specific  heat  and 
vaporization  temperature  for 
target  material 

erg/gm 

DE 

Energy  radiated  away  from  cell 
during  time,  4  t 

erg 

DELP 

(Pi+1  “  Pi)/Pi*  used  in 
forward/backward  differencing 

- 

DELT 

Kronecker  delta 

- 

DELTX 

At/  Ax 

sec/cm 

DELX 

Ax 

cm 

DELX3 

(Ax)  3 

cm  ^ 

DIFC 

Thermal  conductivity 

cal/cm/sec/K 

DIMP 

Impulse  over  the  jth  radial  cell 

dyne-sec 

DK 

DIFC 

DLT 

At 

sec 

DU 

Artificial  diffusion  coefficient 

- 

E(CRNT) 

Mass  specific  energy 

erg/gm 

E (LUX) 

Constant  used  in  finding  DE 

cm/sec 

EDEN 

Electron  number  density 

l/cm^ 

Ell 

Ionization  potential  of  air 

erg 

EI2 

Ionization  potential  of  target 
vapor 

erg 

EL 

Latent  heat  of  vaporization  of 
target  material 

erg/gm 

EM 

Maximum  mass  specific  energy  in 

erg/gm 

Variable 


Definition 


Units 


EMI 

Molecular  mass  of  air 

gm 

EM2 

Molecular  mass  of  target  vapor 

gm 

EN( HYDRO) 

Initial  energy  density/pressure 
in  mesh 

erg/cm2 

EN  ( HOT ) 

Pressure 

erg/cm^ 

ENDEP 

Energy  deposited  in  mesh 

erg 

ENFAC 

Factor  used  to  determine  whether 

EN (HYDRO)  is  pressure  or  energy 
density 

' 

ENV 

Ratio  of  laser  flux  at  radius,  r, 
to  laser  flux  at  center  of  beam 

- 

ERG 

p  v2/2 

erg/cm^ 

F 

Flux  of  hydrodynamic  quantities 
in  cell  n,t 

- 

FACTOR 

Exp  (-ABC  A  X) 

- 

FB 

Laser  flux  incident  on  surface 

erg/cm2/sec 

FBO 

FB/ (EL  +  CTV) 

- 

FNPO 

Flux  of  hydrodynamic  quantities 
in  cel 1  (n  +  1) , t 

- 

FPHMH 

Flux  of  hydrodynamic  quantities 
in  cell  (n-1/2) ,  (t+1/2) 

- 

FPHPH 

Flux  of  hydrodynamic  quantities 
in  cell  (n+1/2),  (t+1/2) 

- 

FWIG 

Flux  of  hydrodynamic  quantities 
in  "zeroth"  cell 

- 

G 

Constant  used  in  finding  electron 
number  density 

- 

GA 

Ratio  of  specific  heats, 7 

- 

GAMMO 

y  -  i 

- 

I 

Index  associated  with  axial 
direction 

- 

Variable 


Definition 


Units 


II 

IM 

IMMO 

IMPO 

J 

JM 

JMMO 

NAT 

NB 

NCYCLE 

NDIF 

NM 

NQ 

NR 

NS 

OMEGA 
P(PRINS) 
P (LUX) 

PC 

PHI 

PHI  2 


Index  defining  spatial  position 
from  front  to  back  of  mesh 

I  defining  maximum  axial  position 

IM  -  1 
IM  +  1 

Index  associated  with  radial 
direction 

J  defining  maximum  radial  position 
JM  -  1 

Variable  which  controls  amount 
of  output 

Index  defining  radius  at  which 
laser  radiation  is  effectively  zero 

Number  of  time  steps  taken 

Variable  which  determines  whether 
energy  density/pressure  is 
differenced 

Index  defining  dimensionality  of  code 

Number  of  hydrodynamic  quantities 
computed 

Variable  defining  radial  position 
in  the  mesh 

Index  defining  cell  position 

EL  +  CTV  erg/gm 

Energy  density/pressure  erg/cm3 

Hydrodynamic  quantities 

Pressure  at  cell  center  erg/cm3 

Factor  used  in  iterative  process 
to  find  temperature 

Factor  used  in  iterative  process 
to  find  temperature 


Variable 

PI 

PINV 

PN 

PNPO 

PNP2 

PNWIG 

PO 

PHMH 

PPHPH 

PRES 

PRESMH 

PRESP 

PRESPH 

PULSE 

R 

R1 

RHO 

RHOI 

RHOO 

RHOOI 

RHOI 

RHO  2 

S 


Definition 

Units 

ir 

- 

l/P 

cm^ /gm 

Hydrodynamic  quantities  in 

cell  n 

- 

Hydrodynamic  quantities  in 
(n+1) 

cell 

- 

Pressure  in  cell  (n+2) 

- 

Hydrodynamic  quantities  in 
"zeroth"  cell 

- 

Ambient  pressure 

erg/ cm3 

Hydrodynamic  quantities 

erg/cm3 

Hydrodynamic  quantities  in 
(n-1/2),  (t+1/2) 

cell 

erg/cm3 

Pressure  in  cell  n 

erg/cm3 

Pressure  in  cell  (n-1/2) , 

(t+1/2) 

erg/cm3 

Pressure  in  cell  n+1 

erg/cm3 

Pressure  in  cell  (n  +  1/2)  , 

(t+1/2) 

erg/cm3 

Total  impulse  over  surface 

dyne-sec 

Radius  from  center  of  beam 

cm 

Root  in  cubic  equation  for 
electron  number  density 

finding 

1/cm3 

Total  density 

gm/cm3 

1/  RHO 

cm3/gm 

Target  density 

gm/cm3 

1 / RHOO 

cm3/gm 

Density  of  air 

gm/cm3 

Density  of  target  vapor 

gm/cm3 

Ratio  of  laser  flux  at  time,  t, 
to  laser  flux  at  time,  zero 
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Variable 

Definition 

Units 

31 

Ratio  of  ionization  rate  to 
recombination  rate  for  air 

- 

S2 

Ratio  of  ionization  rate  to 
recombination  rate  for  target  vapor 

- 

SHAPE 

Same  as  S,  given  above 

SPULSE 

Impulse  over  irradiated  area 

dyne-sec 

SURFP 

Surface  pressure 

erg/cm^ 

T 

Temperature 

erg 

T (GAM) 

Temperature/ 1000 

K 

TEST 

ABC  A  x 

- 

TG 

"Guessed"  temperature/ 1000  for 
iterative  process  to  find 
temperature 

K 

THIRD 

1/3 

- 

THRDX3 

1/27 

- 

TIME 

Time 

sec 

TMID 

Temperature  at  cell  center 

erg 

TN 

Temperature  in  cell  n 

erg 

TNPO 

Temperature  in  cell  n+1 

erg 

TOTE 

Total  energy  in  mesh 

erg 

TP 

Temperature/ 1000 

K 

TPHPH 

Temperature  in  cell  (n+1/2) ,  (t+1/2) 

erg 

TQUIT 

Time  for  calculations  to  stop 

sec 

TWIG 

Temperature  in  "zeroth"  cell 

erg 

V 

Velocity  component 

cm/ sec 

VELSQ 

pv2  /  2 

erg/cm^ 

VO 

Speed  of  sound  in  target  material 

cm/ sec 
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FORMRT < IX *0 16. 10* IX* 14*0 16. 10* IX* I 1> 

RERD  <8  *25I>  RBCM*XULFQ* PULSE  rSPULSE*  TOTE  rENDEP 
FORMRT  <1X*6<016. 10*1X> > 

RERD <8*232>  <FB(I>  *I»1*NB> 

RERD  (8*232>  (CHEW  <I>  *1*1 *NB> 

FORMRT  <1X*7(D16. 10* 1X>> 

DO  63  1*1* IM 
,  DO  62  J*1  *  JM 

NS«NQ*-<l-i+-IM*«J-l>  > 

RERD (8 *23l>  BN <NS*»D  *RN <NS-*2>  *RN <NS+3>  * 
k  BN(NS+4> *RN <NS+3> *BN<NS+6> 

CONTINUE 

CONTINUE 

IF (NCYCLE. EQ . 0>  SOTO  98 
L*NTIME— 1 
DO  61  J*l* JIN 

NS*NQ^(IIN-H*IM-*<  J-l>  > 

T-RN  <NS+G>  / 1 . 380340-13 


UE-AN  <NS+5>  /AN  <NS+t> 

VSQ-O. 

DO  60  NP*l*NW 

VSQ«VSQ+AN(NS+NP>  ++2 
CONTINUE 

VSQ*.  3+VSQ/RN  (NS+1> 

IP(T>  3*3*6 
GAHN0*.3a7 
SOTO  7 

IFCT.LE.3.>  THEN 

GAHHQ*. 4234—.  0768 ♦T 
ELSE 

SAMHC*.  16323+ .  0  0991 7+T 
END  IF 

P*AN  <NS+3>  +<  1-ND  IF>  +6AMH0+ND  IF+<  AN  <NS+3>  -VSQ> 

C*A1  +A2+T+A3+T++2+A4+T++3 

EREF*.  3+<  1000. +C +121. >  +2200.  +R  <  J>  ++2 

TREF*A13+£REF+A1 4+EREF++2+91 3+EREF++3 

A«A3+A6+TREF+A7+TREF+*2+R6+TREF++3 

B*A9+At  0+TREF+R1 1+TREF++2-MH2+TREF++3 

RHOREF*3 . 484D-7+P/TREF* .0 Of 

QDOT  <  J*L>  *<  <A>  ♦+.  6>  +<  <B>  ♦+-.  t>  +<C-.  12l> 

QDDT  <J*L>*t.  0ID+3+DSQRT <RHOREF>+QDQT  <J*L> 
CONTINUE 
DQ  38  J«JIN+t*JN 
QDOT<J*L>«0. 

CONTINUE 
GOTO  99 
NT  IME*NTIHE— 9- 
DD  50  J*1*JN 

WRITE  <9*+>  R(J>  *  (QDOT  <J*L>  *I*1*NTIHE> 
CONTINUE 
CLOSE  <75 
CLQSE<8> 

CLOSE  <9> 

CLOSE  <10> 
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An  attempt  was  made  to  analytically  duplicate  a  set  of 
experimental  data  in  which  enhanced  thermal  coupling  was 
observed.  Analytical  relations  were  developed  to  compute 
the  heat  flux,  into  the  surface  of  an  Aluminum  species,  due 
to  convection  heat  transfer,  target  vapor  condensation,  and 
plasma  reradiation.  From  these  relations,  the  most 
significant  coupling  mechanism  was  identified  to  be  vapor 
condensation.  Convection  heat  transfer  was  found  to  have  a 
less  significant  effect  on  enhanced  coupling,  and  plasma 
reradiation  was  found  to  be  negligible.  The  calculations 
were  performed  by  a  three-dimensional  Lagrangian 
hydrocode.  The  results  of  this  work  are  somewhat 
inconclusive  due  to  the  extensive  computer  resources  used 
by  the  hydrocode;  however,  the  relations  presented  here  do 
seem  to  give  promising  results.  The  work  should  be 
continued  to  investigate  the  effects  of  laser  and  target 
parameters  on  enhanced  thermal  coupling. 
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